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Abstract 


Time-dependent quantum mechanical (TDQM) wave packet (WP) methods have been 
applied successfully in recent years in describing a wide variety of physico-chemical processes. 
In the present study we have used such methods to investigate the transition state (TS) 
resonances and dynamics of elementary gas phase ion-molecule reactions.. 

In chapter 1 of the thesis, we present an overview of the chronological development of the 
TDQM methodology. Its application to a wide variety of problems with a special emphasis 
on TS spectroscopy and reactive scattering is reviewed. 

In chapter 2, the general grid method of solving the time-dependent Schrodinger equation 
(TDSE) is presented. We describe the general set up of a grid and the commonly used 
numerical methods of propagation (in space as well as in time) of the initial WP represented on 
such a grid are reviewed. The stability of various numerical time evolution schemes in presence 
of a negative imaginary potential (NIP), used for damping the WP components near the grid 
edges, is analysed. The Fourier grid Hamiltonian (FGH) method for computing the bound 
state eigenvalues and eigenfunctions is presented. The spectral quantization method (SQM) 
is also presented in this context and its utility in quantifying the resonances is discussed. The 
computation of vibrational state-selected and energy resolved reaction probabilities through 
a time-energy mapping of the reactive flux of the WP, state-to-state reaction probabilities 
through a time-energy projection and state-to-state reaction probabilities from the energy 
resolved S'-matrix elements in the Mpller operator formalism is described at length. 

In chapter 3, we investigate the dynamical resonances in collinear (He, Hj), (He, HD~) 
and (He, DH"*") collisions by analysing their respective TS spectrum on the McLaughlin- 
Thompson-Joseph-Sathyamurthy (MTJS) potential energy surface (PES). The spectra are 
computed by the SQM, through Fourier transformation of the temporal autocorrelation, 
C(t), of the initial WP. The nature of some of the resonances is examined by calculating 
their eigenfunctions and lifetimes. The relation to the corresponding resonant periodic orbits 
(RPOs) in classical mechanics at those energies for the collinear (He, Hj) resonances is 
also established. The vibrational state-selected and energy resolved reaction probabilities 
{P^{E)) for collinear (He, HD"^) and (He, DH"^) collisions are also calculated on the same 
PES by computing the reactive flux in the product channel. The F^{E) values for (He, HD~) 
collisions show a characteristic staircaise-like structure which can be related to threshold 
resonances. The F^{E) values for (He, DH"^) on the other hand are highly oscillatory, in 
keeping with a densely packed TS spectrum. 
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In chapter 4, we examine the possibility of existence of reactive scattering resonances 
in collinear (H~, Ha) collisions by computing the TS spectrum on the Starck and Meyer 
(SM) PES. The nature of the H3 resonances is investigated by computing the corresponding 
eigenfunctions and their lifetimes. We also report the FfiE) values for the collinear reaction 
H~ + H2 (u = 0) -H- Ha + H, on the same PES to confirm the existence' of the resonances. 
The different dynamical behaviors of Hj and H3 are discussed. 

In chapter 5, the TS spectrum of collinear (He, H2) collisions obtained in chapter 3 
is analysed to extract the signature of quantum chaos. Examination of short- and long- 
range correlations in the eigenvalue spectra through a study of the nearest neighbor spacing 
distribution (NNSD), P(s), and the spectral rigidity, A3(X), reveals signatures of quantum 
chaotic behavior. Analysis in the time domain is carried out by computing the survival 
probability, << P{t) >>, averaged over initial states and Hamiltonian. All these indicators 
show beha'vior intermediate between regular and chaotic. A quantitative comparison of << 
P{t) >> with the results of random matrix theory provides an estimate of the fraction of 
phase space exhibiting chaotic behavior, in reasonable agreement -with the classical dynamics. 
We also analyse the dynamical evolution of coherent Gaussian wave packets (GWP) located 
initially in differennt regions of pheise space and compute the survival probability, po-wer 
spectrum and the volume of the phase space over which the WP spreads and illustrate the 
different behaviors. 

In chapter 6, we have presented the results of our preliminary investigation of the TS 
resonances in three dimensional (He,Hf ) collisions obtedned through the spectral quantization 
method. The computer code has been developed in mass-scaled Jacobi (iJ, r, 7) coordinates 
where the spatial evolution of the WP along ( Jh r) is carried out by a two dimensional EFT 
and the evolution along 7 is carried out by a DVR constructed using the Gauss-Legendre 
quadrature. The temporal evolution of the WP is ceirried out by the SO method. 

In chapter 7, a summary of our findings and the conclusions are presented. 

In appendix 1, we have provided some necessary details of chapter 2. 

In appendix 2, energy loss spectra of forward scattered He'*' ions in collision with HCl are 
reported for collision energies (Eiab) in the range, 20-60 eV. There is a -variation of six orders 
of magnitude in intensity for the amount of energy transferred, up to AE ~ 2 eV. Although 
there is very little difference in the extent of AE, for different Eiab, there are noticeable 
differences in the spetra at large AE, indicating propensity for large Au transitions at high 
Eiab. Three dimensional quasiclassical trajectory (QCT) calculations treating HCl as a rigid 
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rotor aad representing the He'^'-HCl interaction as of ion-dipole(quadrupole) in the long-range 
and corresponding to the isoelectronic He'^-Ar in the short-range, could account for much of 
the energy transfer at low Etrans- At higher Etransi rigid rotor model is clearly inadequate 
and slightly better agreement between experiment and theory results when HCl is treated as 
a vibrating rotor. 
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Chapter 1 
Introduction 


Understanding the microscopic phenomena underlying physico-chemical processes is 
the goal of molecular dynamical research. From a theoretician’s point of view such a 
study requires the knowledge of inter- and intramolecular interaction potentials. These 
are either made up by using suitable algebraic functions to model the known interaction or 
evaluated ab initio by solving the electronic Schrodinger equation for the system for fixed 
nuclear configurations under the conditions of Born-Oppenheimer approximation. Once 
the potential energy surface (PES) [1-6] is obtciined the next step is the investigation of the 
internuclear motion on the PES by solving the corresponding equations of motion (either 
classical or quantvun mechanical) [7-16] 

Both experimental and theoretical research in molecular reaction dynamics has xm- 
dergone a revolution during the past few decades. The driving force behind this is the 
technological developments and the emergence of new idecis. Revolution in the experi- 
mental research is largely due to two major technological developments- molecular beams 
and ultrashort (10”^® s = fs) laser pulses. Measurement of steric effect [17] in a reactive 
or non-reactive collision by aligning/orienting [18] the molecules either by using strong 
electric fields in molecular beams [19, 20] or by using lasers in a bulb experiment [21, 22] 
is a great achievement. Rotational cooling of polar diatomic molecules by supersonic jet 
expansion [23] followed by the application of electric or magnetic field perpendicular to the 
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beam flow has created the hybrid, pendular states of molecules [24-30] . Time resolution 
in the measurement up to many orders of magnitude less than the fastest instrumental 
response time has been achieved by losing two synchronised laser pulses- one for excitation 
and another as a probe [31-34]. It has become possible to "clock" the transition state (TS) 
of a chemical reaction (typical lifetime of ~ 10 - 100 fe) in real time and subsequently 
follow the nuclear motion [35-43] . The term transition state above refers to all conflgu- 
rations the system acquires while passing from the reactant to the product VciUey. This is 
slightly different from the traditional definition of the transition state which identifies it 
with the saddle point along the reaction path between the reactants and the products. De- 
velopment of negative ion photodetachment spectroscopy has created an alternative route 
to viewing the transition state directly [44-47]. The neutral unbound transition state is 
obtained in this technique by photodetaching the corresponding bound negative ion. The 
translational energies of the emitted electrons have imprinted on them the ro-vibrational 
structures of the neutral transition state. This has been a long held dream of chemists to 
observe the transition state directly (for a recent review on transition state spectroscopy 
see Ref. [48]). In a recent article [49] Porter commented, "chemists are near the end of 
the race against time". He had initiated the "race" in 1949 through his millisecond flash 
photolysis experiments [50, 51]. 

Much of the early insights into the chemical reaction dynamics had come from quasi- 
classical trajectory studies [52, 53]. Subsequently, the semiclassical methods [54] were de- 
veloped to incorporate quantum features in the classical framework. The time-independent 
quantum mechanical (TIQM) theory [8] of molecular encounters also dominated the field 
until recently. The last few years have seen an intense activity in a time-dependent quan- 
tum mechanical (TDQM) approach [55] to the problem. 

Historically the concept of wave packet (WP) was introduced by Erwin Schrodinger 
(1926) [56]. The thrust behind this was to have a complementary picture of classical tra- 
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jectory in the quantum domain. The correspondence between the two was culminated in 
Ehrenfest’s theorem- " in the classical limit the quantum mechanical expectation values 
behave classically" [57]. The classical limit of course refers to the limi t of zero imcertainty. 
There was no real application of the WP method until Mazur and Rubin [58] solved the 
time-dependent Schrodinger equation (TDSE) for the first time for a model collinear ex- 
change reaction 

A + BC ^ AB -t- C (Rl) 

in 1959. Several years later (1967) numerical solution of the WP motion in one dimension 
was carried out by Goldberg et al [59] and then McCullough and Wyatt [60] studied the 
collinear 

H + H2 H2 + H (R2) 

reaction. They employed an implicit (Crank-Nicholson) scheme for the time evolution of 
the WP and a finite difference scheme for its spatial evolution. Their method involved 
tedious matrix inversions and was fovmd to be computationally intractable in higher di- 
mension. Later on (1978), Askar and Cakmak [61] proposed an explicit time propagation 
scheme, which turned out [62] to be a re-discovery of the earlier work of Harmuth [63]. 
Their method was devoid of the tedious matrix inversion. Feit et al [64] in 1976, had 
introduced n Fourier based pseudo-spectral scheme to study the propagation of the laser 
beam through atmosphere. They proved the appKcability of their method to solving the 
TDSE numerically on a grid to compute the energy level spectrum of the Henon-Heiles 
system [65] . The Fomier method offered an accurate way to evaluate the action of ki- 
netic energy operator on the WP by switching over to momentum space. This wets the 
first breakthrough in the numerical technology of solving the TDSE. The accuracy of the 
method resulted in abandoning the finite difference scheme which was not very stable. Feit 
et al [65] used a split-operator scheme for the time propagation. Subsequently, Kosloflf and 
Kosloff [66, 67] reiterated the usefulness of the fast Fourier algorithm combined with a 
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second order differencing scheme (explicit scheme of Harmuth) to solve the TDSE. Light 
et al [68] developed a general quadrature based pseudo-spectral scheme (generalised dis- 
crete variable representation (DVR)) to carry out the spatial propagation of the WP. It 
was found to be more promising to evaluate the angular kinetic energy part where the 
usual FFT method faces a singularity at 6=0 and -tt. Tal-Ezer and Kosloff [69] developed 
a Chebyshev polynomial scheme for the time evolution of the WP and found it to be 
the most accurate. The Chebyshev scheme combined with the FFT and/or DVR method 
is often the method of choice for solving the TDSE, particularly, for systems exhibiting 
long-time dynamical behavior. Bisseling and Kosloff [70] showed the applicability of the 
Hankel transform in the place of the FFT and pointed out its usefulness in dealing with 
higher dimensions. The short iterative Lanczos scheme was introduced [71-74] as another 
alternative for the time evolution and found to be quite useful in dealing with the explicitly 
time-dependent Hamiltonians. A comparative study of different numerical time evolution 
schemes with respect to their accuracy and stability in solving the TDSE has been carried 
out by various researchers in the literature [75, 76]. A brief summary of various numerical 
schemes for solving TDSE is presented in chapter2. 

The TDQM approach provides a realistic view of chemical reactions in a theoretical 
calculation. It allows .one to monitor the progress of the reaction event in time by looking 
at the snapshots of the WP on the potential grid, for example. It has a high degree*-of 
predictive value, and it also suggests means of controlling chemical reactions [77-86]. Prom 
a practical computational point of view, the TDQM approach offers yet another advantage 
in that, a single wave packet propagation yields the desired result like reaction probability 
at a range of energy values. TIQM calculations, in contrast, pose a serious problem as a 
high resolution in energy requires the use of a dense grid in the energy scale, which results 
in a substantial increase in computational time. Finally, near the onset of the collision- 
induced dissociation one runs into the basis set problem in the TIQM calculations, the size 
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of the basis set becoming very large to take care of the contribution from the dissociation 
continuum, whereas, TDQM methods pose no such problem. 

During the past few years the TDQM techniques have been used to solve a variety 
of problems, half-collisions as well as full collisions. In a half-collision study one starts 
from the transient and study its decay into products. The classic example of this is 
photodissociation and predissociation processes [87-96] In photodissociation, the molecule 
is prepared in the specified rovibrational state of the ground electronic state and then 
taken to the excited repulsive electronic state, which may be coupled with other excited 
electronic states, through the action of the transition dipole matrix elements and from 
where it fragments into products. The quantities of interest here are the photodissociation 
cross sections (partial and total), anguleir distribution of the photofragments, and the 
branching ratios if there is more than one product channel. In predissociation, the species 
is excited from a bound (ground) to a bound (excited) state and before the excited species 
emits light and returns back to the ground state it dissociates due to the coupling of the 
upper bound state with a repulsive state. Recently the idea has been extended to compute 
transition resonances in chemical reaction [97-102] . Resonances in a chemical reaction 
arises from the presence of metastable states in the transition state region of the PES. 
Characterization of resonances and their influence on the dynamics have played a major role 
in research in the area of reaction dynamics [103-106] . Some of these resonances have been 
identifled in the photodetachment experiments of Neumark and coworkers [44-46, 107-110]. 
Their experiments combined with the coupled channel h 3 rperspherical reactive scattering 
calculations of Schatz and others [111-116] have been utilised extensively in understanding 
the topology of the PES for these systems. In addition the TDQM approach has been used 
in several other applications, such as multiphoton dissociation processes [117-119], mode 
selective decay [120-131], vibrational predissociation of van der Waals complexes [132, 
133], Raman scattering [134-136], computation of the cumulative reaction probabilities 
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[137, 138], computation of the reaction probabilities by computing the scattering matrix 
in the M0ller operator formalism [139-143] , investigation of the nature and dynamics of 
the WP prepared by an ultrashort laser pulse [119, 144-153] , etc. 

In a full collision study the reactants are prepared selectively in the desired internal state 
and then allowed to interact to lead finally to products. In addition to a study of reactive 
scattering processes [154-167] , inelastic [168] and gas-surface scattering [169-186] processes 
have also been treated by the TDQM method. The methodology has been extended to 
four atom systems [166, 187-190] recently. Several excellent reviews and monographs [191- 
199] on the TDQM approach to reaction dynamics have appeared in recent years. Some 
aspects of the TDQM methodology are presented in details in chapter 2. We have used the 
TDQM method to investigate transition state resonances in collinear (He,Hj), (He, HD-^) 
and (He, DH"*") collisions. The results are given in chapter 3. For the latter two systems, 
reaction probability values as a function of relative translational energy of reactants are 
also presented in chapter 3. Results of the investigation of TS resonances in collinear 
(H“, Ha) collisions form the subject of chapter 4. Results of chapter 3 are analysed from 
quantum chaos point of view, in chapter 5. The methodology and preliminary results of 
a three dimensional investigation of (He, H^) collisions are given in chapter 6. Finally a 
summary of our findings and conclusions axe presented in chapter 7. 



Chapter 2 
Methodology 


2.1 Introduction 

In time-dependent quantum mechanics the wave function, ^ at time, t, is evaluated by 
solving the time-dependent Schrodinger equation (TDSE): 

ih— = m (2.1) 

where I’d- is the Hamiltonian operator. In order to solve this equation numerically 
it is necessary to represent the initial state of the system and follow its evolution in time. 
The first step in this endeavovu: is accomplished by constructing a discrete Hilbert space 
in the form of a grid. This is often referred to as the grid method in the literature. 

The molecular encounter takes place in the position-momentum and the time-energy 
space. The discrete grid is usually constructed in the coordinate space, with each grid point 
being characterized by a finite value of the interaction potential. Once this grid is set up 
the overall solution follows the next three steps: (i) preparation of the initial wave packet 
(WP) = 0)) (depends on the system and the desired observable) and representation 
of this WP on the grid, (ii) its evolution in space and time and (iii) final state analysis. In 
this chapter we will briefly review various numerical techniques for performing steps (ii) 
and (iii). Step (i) will be described while dealing with a particular system. 

One serious problem in the TDQM calculation is that as time progresses the WP 
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gradually reaches the grid edges and undergoes spurious reflections resulting in interference 
between the outgoing and the reflected components. One way to circumvent this problem is 
to use a very large grid which would delay the WP reaching the boundaries. But this results 
in a tremendous increase in the computer memory and time requirements particularly in 
higher dimensions and therefore is not feasible in most practical applications. In the late 
eighties Neuhauser et al [200, 201] have introduced a negative imaginary potential (NIP) 
of the form — iVb, in addition to the real potential of the system at the last few points 
nem the edges of the grid. The use of NCPs is well known in optics. They provide a 
convenient way to damp out the WP components nem the grid boundary and thereby 
prevent the unphysical reflection. This alleviates the necessity of very Imge grids and 
minimizes the problem of computer memory and time. Although NIPs are often used in 
TDQM calculations, their use causes problems in some of the numerical time evolution 
schemes as the Hamiltonian becomes non-Hermitian in their presence. This, naturally 
raises questions about their use in time-dependent quantum moieculeu scattering. We 
illustrate this point at length later in the present chapter. An alternative and a better way 
to prevent the WP reflection at the edges by using real damping functions is also discussed 
in this chapter. 

2.2 The Spatial Grid 

Let us asstune that the coordinate space is discretized into a set of N discrete points. If the 
spacing between two successive points is Ax, then the eigenvalue of the position operator 
X at each point on the grid is given by 


Xi = {i — l)Ax, 


i = l,...,iV 


( 2 . 2 ) 
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The corresponding eigenvectors are denoted by [xj >. The orthogonality and the com- 
pleteness relations on this discrete grid are given by: 


2 Air(xj|l,) = Sif 


(2.3) 




and 


N 


^ |a:i > Ax < Xj] (2.4) 

i=i 

where is the identity operator. The functions represented at the grid points are given 
by 

{xi\(f>) = o(xi) (2.5) 

The continuous normalization integral, J^(p*{x)o{x)dx = 1, on this grid transforms to a 
discrete sum: 

N 

Y2(f>*{xi)<j){xi)Ax = 1 (2.6) 

i=l 

The maximum length of this grid along the spatial coordinate x is, T = NAx. This length 
determines the spacing between two successive points in the reciprocal k space: 

27T 


Ak 


NAx 


(2.7) 


In the k space the grid is centered at zero and ail other points are distributed symmetrically 
on either side. If Pmax{— represents the maximum momentum in the k space then 

the N points are distributed in the interval {—pmax, — , 0, ....iPmax}- Hence the total length 
of the grid in this space is 2|pmoi|- Now Ax can be written as: 


Ax = 



( 2 . 8 ) 


Therefore, the total volume represented in the phase space 


Volume = 2Lpmax = Nh 


(2.9) 
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is decided by the number of grid points. The maximum energy represented on the grid 
through this discretization is given by: 


E = T 4-F 

^m€ix -‘•max ‘ ^max 

rmca , xr 

“ 2m ^ 




+ K. 


2m ' 

= Y — + y 

/- ^ / A \9 ' max 


( 2 . 10 ) 


^ 2m(Aa:)2 

where Ymax is the maximum value of the potential energy represented on the grid. The 
minimum energy £'m»n(= Kiin) is equal to the minimum value of the potential energy. Our 
subsequent discussion will follow upon this general set up of the discrete grid. 


2.3 The Spatial Evolution of the Wave Packet 


In order to evaluate the action of JT on ^ we need to evaluate separately the action of the 
kinetic (T) and the potential (T^ energy operators on it. Fis local in the coordinate space, 
so its action on ^ is only a multiplication of its magnitude with the value of ^ at the grid 
points: 

T^a:)^(xi) = (2-11) 


The kinetic energy operator {T = ^ = ^ non-local in the coordinate space and 

the evaluation of constitutes a key step in quantum molecular dynamics. This is done 
either by a semi-local approximation in the finite difference (FD) scheme [60] or in the 
local representation by a suitable collocation technique [202]. 

In the FD scheme is evaluated on the coordinate grid, for example, through: 


TpD'^iXj) 


2m 

[ ^(a:;+i) + ^(xj-i)-2^(Xj) 

2m 


( 2 . 12 ) 
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This is referred to as the 3-poin.t FD method. Here the wave function is approximated 
locally by a polynomial which results in a power law convergence of the method, (Aa;)'‘. 
The semi-local description, unfortunately, does not conform to the commutation relations 
of quantum mechanics. The error accumulates exponentially on repetitive application of 
Ton 


2.3.1 The Collocation Technique 

The collocation technique is a global optimization technique pioneered by Gauss [203]. 
The basic idea imderiying it utilizes the concept of discrete Hilbert space. Continuous 
functions are approximated at each point in this space by a discrete sum in terms of 
a finite functional basis set. The basis functions (ffn(xi)) at various grid points Xj axe 
connected through appropriate expansion coefficients (on) to approximate the continuous 
fu n ction at the grid points: 

.v-i 

^(Xi) = f (Xi) = ^ an^n(Xi) (2.13) 

n=0 

where N is the size of the basis set. This is also known as the pseudospectral approximation. 
As described by Kosloff [202, 204, 198], the functional space in this method is viewed 
through a set of portholes, which is restricted to a local picture, and if the portholes are 
dense enough the full functional space is mapped on to the discrete space by interpolating 
the observed data by utilizing the continuous properties of the functional space. In matrix 
notation the above equation can be written as: 


^ = Ga (2.14) 

This is a set of coupled linear equations with = ^(xi) and G„i = gn(sci). The solution 
for the expansion coefficients for a set of linearly independent y„(xi) becomes 


a = 


(2.15) 
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This is the general collocation scheme which provides a global mapping technique provided 
by the functional basis connected through the expansion coeflScients, On, to the spatial grid. 
The scaleir product of two functions on the functional space is given by 


(®(X)|4(X)> = 

nm 

where Snm is the overlap matrix given by 

(2.16) 

= / dxg*„(x)g^{x) 

J D 

with 

(2.17) 

^(®) = J2‘^9nix) 

n 

and 

(2.18) 

^(a:) = II bmgmix) 

m 

(2.19) 

The collocation scheme (vide supra) becomes simplified when the functional basis is chosen 

firom an orthonormal set. In this case 

N-l 

H 9ni^i)9nixj) Sij 

n=0 

Now left multiplying Eq. (2.13) by ^^(si) and then summing over i results in: 

(2.20) 

On = J^^(Xi)g*^{Xi) 
i 

(2.21) 

which shows that the coefficients a„ are the discrete functional transform of 

condition Eq. (2.17) becomes: 

Under this 

•Snm ~ bnm 

(2.22) 


and Eq. (2.16) 


n 

= f; 

1=0 


(2.23) 
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The matrix G becomes unitary, and the Hilbert space becomes a discrete vector space 
with a unitary transformation between the discrete collocation points Xi and the discrete 
functional base a„. A special case of this orthogonal collocation scheme is the Fourier 
method. 


2.3.2 The Fourier Method 

In the Fourier method the orthogonal functions {^A:(aj)} are chosen as the plane wave basis 
functions [64, 202, 204, 198]: 

gk{x) = k = -{N/2 - 1), 0, ...A/2 (2.24) 


Expansion of a wave function in this basis is given by: 


iV/2 

h=-(N/2-l) 


(2.25) 


which on inversion using the set of equidistant discrete sampling points { 0 :^} gives the 
discrete Fourier expansion coefficients {at}, the amplitude of the wave function in the 
reciprocal k space: 

= (2.26) 

The factor ^ arises from the following orthogonality relations satisfied by the Fourier basis 
set: 


E 9k(o:i)9; {zi) = SkiN, \k - i\ < N (2.27) 

i=l 

Eqs. (2.25) and (2.26) constitute the discrete versions of the continuous Fourier transfor- 
mation 

^x) = -4= r ^{k)e^^^dk (2.28) 

V 27r J —00 

and 


1 f°° 


^k) = / ^x)e-^'^dx 

VStT J-oo 


(2.29) 
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Operationally, one can switch back and forth between the two reciprocal discrete Hilbert 
spaces (position and momentum or the time and energy) through Fourier transformation. 
It has been shown [66] that in doing so the Hermitian properties of the operators and their 
commutation relations are preserved. The domain in which the transformation takes place 
has periodic boundary conditions. 

The kinetic energy operator f (=^) is local in the momentum space and hence is 
diagonal in k. Therefore the 2^(x) operation is carried out by Fourier transforming ( 
FT"^) ^(cc) to followed by a multiplication by the kinetic energy spectral value at 

the discrete momentum grid points 

- ft2.2 _ t2T2 

IW(x)=— = — a. (2.30) 

and then bringing it back to the coordinate space by an inverse Fourier transform (FT“). 
The domain in which the transformation takes place requires the periodic boundary con- 
dition. The evaluation of the derivatives becomes exact for wave functions which obey this 
periodic boundary condition. The quantum mechanical commutation relations remain pre- 
served for periodic potentials and functions under these potentials remain localized in the 
phase space box and the amplitude of the function becomes zero at the boundciry of the 
box. For wave packets this boundary condition is met when they are band-limited func- 
tions with finite support [202]. The support of a function is defined by a set for which it has 
non-zero value. The wave functions are not usually band-Hmited with finite support since 
they can not be confined in both coordinate and momentum spaces. Wave packets which 
are semilocalised wave functions satisfy this requirement. The use of absorbing bound- 
aries also forces wave packets to meet this requirement by reducing their amplitude to 
zero at the boundaries. Another important aspect which has made the Fourier transform 
an attractive tool in molecular dynamics is the availability of the fast Fourier transform 
(FFT) algorithm resulting firom the work of Cooley and Tukey in the mid ’60s [205] and 
others [206, 207] and the fact that it could be adapted to the vector and parallel computer 
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architectures in more recent times [208] The FFT algorithm scales with N as 0{NlogN). 

It is shown by Kosloff [66] that the numerical dispersion in the FFT scheme is much 
less than that in the case of the FD scheme. The kinetic energy operator in the latter 
scheme, on Fourier transformation to the A:-space gives [202] 

2(cos(A; A x) — 1) 


TMk) 


2m 


(Axy 


2sin{kAx/2) 
2m Ax 


t2 


(2.31) 


The FD spectrum obtained from the use of the above equation deviates considerably 
from the Fourier spectrum as k increases. The deviation is shown [202] to be comparatively 
less in the case of higher order FD spectra. Kosloff [202] argues that acceptable accuracy 
can be obtained by using at least ten points per wave period in the FD method. This means 
that one hzis to use ten points per unit volume in phase space. Moreover, as mentioned 
above, the semi-local approximation in the FD scheme does not preserve the commutation 
relations of quantum mechanics. The FD algorithm scales as 0{N^) which is inefficient 
for higher dimensional problems. All these factors established the superiority of the FFT 
scheme over the FD scheme. 


2.3.3 The Hankcl Transform 

The fast Hankel transformation (FHT) algorithm was outlined by Siegman [209] and also 
by Talman [210]. Bisseling and Kosloff [70] introduced the FHT in molecular dynamics in 
the context of treating the radial part of the Laplacian operator in polar coordinates. The 
Hankel transformation of order n of a function F, from the r space to the corresponding 
reciprocal k space in the interval [0, oo] is defined by: 

F{k) = rF{r)J^{kr)dr, k>0 (2.32) 

where J„’s are the order Bessel functions of the first kind. This is also recognised in the 
literatmre as the Foxurier-Bessel transformation. Operationally, FHT is performed through 
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a logarithmic change of variables [209, 210], which substitutes r = roe ^ and k — in 
the above equation and results in: 


F{kQd^) =rl r° e ^^F[roe ^)Jn(rokoe^ ^)dy 

J — OO 


(2.33) 


In the discrete Hilbert space F is sampled at N (chosen in powers of two) uniformly spaced 
(A) discrete grid points (x = = mA;i^Tn = 0, 1) and the above integral is 

evaluated, as a discrete sum: 


Fikoe^^'-) ^ rlA ^ e-^-^i^troe-”*^) Jn(roA:oe('-"*)^) (2.34) 

m=0 

The right hand side of the above equation represents a convolution of two fimctions, which 
is carried out using the convolution theorem: "the Fourier transform of the convolution of 
two functions is equal to the product of their Fourier transforms" [211]. The evaluation of 
the continuous (non-circular) convolution by the discrete (circular) Fourier transformation 
is done by padding the values of F with N zeroes by setting: 


bm = 


0 for m = —N, , — 1 

foj. ^ = Q, N- 1 


and 


Q = A for i = — iV, ,N— 1 


(2.35) 


(2.36) 


This changes Eq. (2.34) into a 2N- term discrete circular convolution 


ai=T,^ZlbmCi.m fori = -iV; ,N-1 


(2.37) 


The 2iV’-point FFTs are separately performed on b and c and the results are multiplied and 
an inverse FFT is perfonned on the product. The first N values of the result are ignored 
and the remaining N veilues of = F{kod^) for i = 0, ..., iV — 1 are the approximated 
values of the Hankel transform oi F. 

Both real and complex functions can be evaluated by this method. The parameter 
ro = Tmax Is the maximum value of r and k) = kmin is the minimum value of k. The other 
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parameters are chosen such that the domains of and k{kmin^ kmax) contain the 

function F. The sampling points in the discrete space are chosen such that one point falls 
in between two zeroes of the Bessel function which occur at a separation of tt. This 
leads to: 

lroA:o(e(^“^)^ - e(^-2)^)| < tt (2.38) 

or approximately 

< tt (2.39) 

or 

^maxkjnax^^ ^ (2.40) 

where kmax = A:o(e^^“^^^). The above restriction provides a criterion for choosing A. It is 
suggested by Bisseling and Kosloff [70] that the distance between two sampling points of 
F{— e“”*^ — should be less than the minimum distance between two successive 

zeroes of F{= Tr/kmax) represented on the logarithmic grid. This ensures a sufficiently 
dense sampling for F. 

The operator [7 ^ is evaluated by performing FHT of order n followed by 

multiplying the result by —k"^ and again performing the order FHT. The FHT is its 
own inverse. An FHT of orde zero is performed for operators which do not contain terms 
like 

The FHT method is prescribed for systems with a radially sy mm etric potential. For 
such systems the dimensionality of the problem can be reduced by switching over to the 
polar coordinates. The extra computation involved in the irrelevant portion of the grid 
(which has no value from the dynzimical point of view) in Cartesian coordinates is reduced 
in polar coordinates. The method has been applied by Bisseling and Kosloff [70] to compute 
the energy eigenvalues of 2D and 3D harmonic oscillators and also for studying H + 
H2 reactive collisions. Their estimates show that both FFT and FHT methods are of 
comparable accuracy. However, the CPU time is estimated to be less in case of the FFT 
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and it is expected that FHT will be faster for systems with a small skewing angle and 
where the coupling between the radial and the angular part is weak. It must be added 
that the FHT algorithm scales with N as NlogN. 

2.3.4 The Discrete Variable Representation 

The discrete variable representation (DVR) is another example of the orthogonal collo- 
cation method which uses specific functions and points. The method uses an orthogonal 
transformation between the DVR and the finite basis representation (FBR) and vice versa. 
In this transformation the Hermitian property of the operators is preserved. The DVR is 
basis localized at the discrete eigenvalues of the coordinate operator, whereas the FBR is 
the basis consisting of N square-integrable functions of the continuous coordinates. The 
latter is also known as the spectral representation. The method which utilizes a different 
representation of the wave function is called the pseudospectral method. Both the Fourier 
and the DVR are methods based on the evaluation of operators in their respective local 
representations. The DVR was originally developed by Harris et al [212] to evaluate the 
potential matrix elements in the harmonic oscillator basis (the Gauss-Hermite points) by 
using the DVR-FBR transformation. Subsequently, Dickinson and Certain [213] showed 
that classical orthogonal polynomial basis functions, (e.g. Hermite polynomials) corre- 
spond to Gaussian quadrature and an orthogonal transformation exists between the two 
representations; N quadrature points and iVbasis functions. Light and coworkers [68] devel- 
oped the DVR in the context of scattering studies. Independently, Shizgal and Bladcmore 
[214] developed DVR under the name of discrete ordinate method. Although initial inter- 
est in the DVR was to describe the internal degrees of freedom in the scattering problems, 
it has been extended to a wide variety of other problems now. The method is also found 
to be quite advantageous for multidimensional problems. 

As was pointed out earlier, the DVR utilizes the idea behind the pseudospectral or 
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collocation method of using a mixed representation of the function ^(x). If a grid repre- 
sentation of {^(xi),i = is used then the action of any local operator (j4) 

on it is evaluated by multiplying its value j4(xj) at the grid point with the corresponding 
value of the function ^'(xi) at that point. On the other hand, the action of a non-local 
operator (B) is evaluated in the finite basis set representation {<?in(a:)} where the operator 
B is local by switching from the grid based representation to the basis set representation. 
It is shown by Dickinson and Certain [213] that an isomorphic transformation exists when 
the transformation matrix is unitary. This is achieved by using the orthonormal (L^ ) basis 
set of some classical polynomials for {^„(Xa)} defined on a set of discrete quadrature points 
(xa, a = 1, ..., N) with the associated weights Wa- One can, in this discrete representation, 
write: 

(^i^nka) = ^n(Xa) (2.41) 

In pseudospectral approximation, ^(x) can be represented as in Eq. (2.13) through [215]: 

'^(^) (2.42) 

n 

The orthonormal property of ?i„(x) is utilized to evaluate c,i through 

Cn= f (j)*^{x)'^{x)dx (2.43) 

Jd 

If the sampling points are chosen as the zeroes -of the ^;^+i(x) polynomials then the theory 
of Gaussian integration enables us to write the orthonormaiity condition of <j)n{x) at the 
quadrature points as 

Y,W„(j)*^{Xo,)(l>m{Xa) = Snm (2-44) 

a 

This is used to write the inverted collocation relation at the N discrete quadrature points 
as 


(2.45) 
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using which. ^ (x) can be evaluated at the quadrature points through 

n a 

= E (2.46) 

a 

where 


Sa(x^ yj'^a. (^a)^n(^) 

(X 

= (2.47) 

n 

"f 

Tni is the transformation matrix from the N dimensional spectral representation (FBR) 
to the N dimensional quadrature representation of the discrete Hilbert space (DVR): 

tJ; = (2.48) 

The orthonormality condition mentioned above in Eq. (2.44) can be rearranged to show 
that T is column orthogonal [(T ^T)^^ = •S'nm = < 571771 ]. In tlie discrete representation the 
overlap matrix is given by [215] 

= (a|/3) 

= / Sl,{x)60{x)dx 

Jd 

n 

= (TTt)„^ (2.49) 

which follows from the orthonormality of the spectral bcisis. Since T is column orthogonal, 
the above overlap matrix A is idempotent [A^ = A] and it will have elements either zero 
or one. Moreover, tr{T ^T) = N, the number of basis fxmctions and A = T will have 
N eigenvalues equal to one, and therefore is an identity matrix. Hence, T,il is orthonormal 
and unitary (T = TT ^ = 1). Therefore T corresponds to the transformation matrix 
from the DVR to the FBR. {5a (a^)} are the spectral basis functions rotated by the unitary 
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transformation matrix dependent on the continuous variable x, and eire the basis functions 
mapped onto the discrete Hilbert space. These are discrete analogs of the Dirac delta 
functions emd hence eire the eigenfunctions of the coordinate operator with eigenvalues 
equal to the quadrature points, Xa 

( q :| x |/ 3 ) = Js*^{x)x 60 (x)dx 
C 

(2.50) 

If the potential matrix is diagonal in the DVR, then its elements at the discrete points 
cire given by 

{a|i/®™|/3) = (2.61) 

The DVR of the Hamiltonian can be constructed as 

nm 

= (T (2.52) 

Whitnell et al [216, 217] constructed the DVR in manydimensions as a direct product 
of one dimensional DVRs by successive diagonalization, truncation and recoupling [218]. 
This becomes particularly simple when an orthogonal coordinate system is used, in which 
case the Hamiltonian does not contain mixed derivatives. 

^ 4“V(xQt, j/^, (2.53) 


where the kinetic energy matrices all have the form 

K = T 


(2.54) 
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In higher d i m ensions the matrix multiplication will have an iV^ scaling relation, where 
N is the dimensionality. This is to be contrasted with the NlogN scaling relation found 
for FFT. Colbert and Miller [219] have shown the NlogN scaling can be achieved in multi- 
dimensional DVR by using Cartesian coordinates without requiring emy faster algorithm. 
Their argument was that in Cartesian cordinates the DVR of the kinetic energy operator is 
sparse, i.e, it can be written as a sum of one dimensional kinetic energy operators. Corey et 
al [215] have pointed out subsequently that even a non-Cartesian coordinate representation 
of the Hamiltonian matrix leads to an NlogN scaling relation. Their argument is that the 
scaling relation is a result of neither the separability of the kinetic energy operator nor 
the sparseness of the Hamiltonian matrix, but due to the fact that any multi-dimensional 
Hamiltonian can be written as a sum of products of functions of the local (coordinate) 
operators and two derivative operators which are nonlocal only in a single coordinate. 
Bramley et al [220] have pointed out that any multidimensional operator is factorizable if 
all its off-diagonal terms couple to a single coordinate. 

A wide variety of classical polynomials, viz., Jacobi, Legendre, Laguerre and Hermite 
polynomials are used as the finite basis set. The grid is defined by the Gaussian quadrature 
points and the discrete inner products in the FBR aire evaluated. It has been pointed out 
by Light and coworkers [68] that the choice of quadrature points which are the zeroes 
of the {N -h l)th term of the orthogonal polynomieil eliminates the contamination of the 
numerical quadrature firom the non-orthogonal components of the projected function which 
lies outside the Hilbert space. They have also defined a unitary transformation between 
the FBR and an arbitrary set of coordinates. There are examples of many other DVRs. 
For example, Muckerman [221] has constructed the DVR basis functions which satisfy the 
Kroenecker delta condition for a general quadrature consisting of a set of iV points and 
weights. He has shown that by employing the Fourier basis set a periodic DVR results and 
leads to the equally sampled points of the Gauss-Chebyshev quadrature. Manolopoulos et 



23 


al [222] used a Lobatto shape quadrature and Lobatto shape functions. They defined a 
DVR basis, which satisfies the Kroenecker delta relation by 

N-i-l 

Qiix) = n {ij = 0, l,j # i) (2.55) 

i=o 

The advantage of this DVR is that the end points of the coordinates can be included 
explicitly with appropriate boundary conditions. Shapiro and coworkers [223] and Clary 
[224] have proposed a DVR that is generated by diagonalizing the coordinate operator. 
However, the basis functions are not polynomials and the DVR does not correspond to a 
Gaussian quadrature. The DVR basis functions also do not satisfy the Kroenecker delta 
relation. Colbert and Miller [219] proposed a novel basis free DVR. It consists of an 
infinite number of equally spgiced points from which one extracts the finite subset desired. 
The kinetic energy operator in this method is defined either in the infinite range and 
approximated by an infinite order finite difference approximation leading to the normal 
Chebyshev DVR or as the finite limit of the particle-in-a-box functions or Chebyshev 
polynomials. Unfortunately, the finite range definition in their scheme faces normalization 
problems. 

The FFT scheme is numerically ineflB.cient in solving the TDSE in the polar coordinate. 
This is because the kinetic energy operator in the polar angle B contains a (1/ sin^ 6) term 
which leads to singularity in the discrete angle space for 0 = 0 and x. One can deal with 
this situation either by performing a discrete cosine transform or by using a DVR. 

The rotational kinetic energy operator is given by [225] 



' 1 d . £ 

sin 9 89^^ 89 



(2.56) 


The discrete cosine transform was proposed by Dixon [226] in which the singularities at 9 
= 0 and X are avoided by sampling 9 at a half integral regular mesh of points between 0 
and X. The Fourier amplitudes F{j) in the angular momentum space are then computed 
by performing a discrete quarter- wave Fourier transform [227]. 
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This situation is dealt with in the DVR by performing the kinetic energy operation 
in the angular momenttun space (FBR). In this case the DVR points are chosen as the 
abscissa points in Gauss-Legendre quadrature, with associated weights Wa. The FBR is 
defined in terms of the normalized associated Legendre polynomials, QjmiXa), of degree 
— Tn|. These functions are orthonormal in the interval [—1,1] and the orthonormality 
condition is given by 

©/m(Xa)©im(Xa) = “ Xa)'”‘'-P/m(Xa)-P;m(Xa)dXa = % (2.57) 

where Xa = cos^^ and (1 — ^rot, is diagonal in this basis and it results in: 

tatejmiXa) = ^ ^ 0;>n (Xa ) (2.58) 

The treatment of the anguleir coordinate in the DVR representation is described in de- 
tail by Leforestier [73, 228] and also by Corey and Lemoine [229]. These authors have 
formulated the scheme utilizing the generalised DVR formulation of Light and coworkers 
[68]. Their method has been adapted by Balint-Kurti and coworkers [230] for studying 
the photodissociation of HOCl and also for computing the state-to-state cross sections for 
3D (Li,FH) reactive collisions [231]. The method uses the same 9 grid for the different m 
components of the wave function. The imitary transformation matrix from the DVR to 
FBR is given by 

= v^0>.(Xa) (2.59) 

The grid wave function is defined as 




( 2 . 60 ) 


and r,.ot^”'(0„) is evaluated at the discrete grid point by 




;=|ml 


(2.61) 
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where in the DVR is mapped onto the FBR through the transformation matrix 

T’J^ and multiplied by the diagonal matrix elements of Trot, , and then transformed 

back to the DVR through the Hennitian conjugate of T^. The number of an gular 

basis functions jmca = {Ng — 1) where Ng is the number of grid points along 9 for each m 
is found to give good performance although different authors have used different numbers 
of angular basis functions. 

2.4 The Time Evolution of the Wave Packet 

2.4.1 The Time Evolution Operator 

The solution of the TDSE is given by 

^(0) (2.62) 

where ^(0) is the initial wave function and ^(t) is the wave function at time t and P is the 
time ordering operator [232]. For an explicitly time-independent Hamiltonian the above 
equation simplifies to; 

^(t) = (2.63) 

The exponential operator in the r.h.s. of the above equation forms a continuous dynamical 
group where time i is a parameter, and is known as the (time) evolution operator denoted 
by tr(^,to). For *0 = 0, 

(2.64) 

In actual computation, t is sliced in smaller steps of length, Ai, and the entire time 
evolution is accomplished through: 

Nt-\ 

n &[(n + 1) A t, n A t] (2.65) 

n=0 

where Nt is the total number of time steps and At = t/Nt. U{t, to) is a linear operator and 
is unitary: 


T 

Pexp -i/h f H{t')dt' 
Jo 


4 * ^ 4 * 

UU* = 17 ' 17= 1 


( 2 . 66 ) 
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The exponential operator is approximated in various ways. For example, by the second 
order differencing (SOD) scheme, the spht operator (SO) method, the Chebyshev polyno- 
mial (CP) scheme or the short iterative Lanczos (SIL) scheme. Each of these schemes is 
briefly reviewed below and special attention will be paid to the stability of these schemes 
in the presence of a negative imaginary potential (NIP). 

2.4.2 The Second Order Differencing (SOD) Scheme 

The time evolution operator and its adjoint for a time step of length {At) can be written 
as; 


Therefore, 


U{At) = 


U^At) = 

(2.67) 

^(f + At) = U{At)<S{t) 

(2.68) 

^'(t- Af) =ut(Ai)^(f) 

(2.69) 


Using the above equations one can write 

^{t -t- At) - '5(f - At) = - e’^'^*‘^^)^(f) 


(2.70) 


Expanding the exponential terms in the above equation in Taylor series and keeping only 
the terms up to first order in At, the equation becomes 

^(f 4- At) = ^(f - At) - ^^^^^{t) (2.71) 

Thus the SOD scheme [61, 66] evaluates ^(f -I- At) from its values at time t and t — At. 
Computationally, to initialize the iteration process, the value of 'S' at the first step is 
computed either by a fourth order Rimge-Kutta scheme or by a higher order Taylor series 



27 


expansion of U{At). The above equation can also be written as: 


h 

= Qm 

The scheme will be unitary if QQ^ = = 1. Prom the above it follows that 

OQt = 


iMt/h SiHAi 

-iHAt/h . ^iHAt 

ft 

h 


= l + (4/6)-^ 


(2.72) 


(2.73) 


neglecting higher order terms in the expansion of the sine function. The conservation of 
the norm can be checked through: 


< ^(t + Af)|^'(t + At) > 


< ^(<)|^(t) > +(2/3)(At)^ < J^'P(t)|.^^(t) >(2.74) 


Therefore it becomes clear that by choosing At sufficiently small the unitarity can be 
maintained throughout the time evolution. It can be shown that [66] the algorithm is stable 
if At << hfEmax ; where Emax{= Tmax + Imax) Is the maximvun total energy represented 
on the grid. In practice, a choice of At < hf^Emax is found to give good results [66]. 

Since the Hamiltonian for the present case commutes with itself, ^ = 0 aind 

energy is conserved throughout the time evolution. 

If a negative imaginary potential (P =— iVo) is added to the real potential of the system, 
the Hamiltonian becomes complex = H+T) and Q' can be written as 

Q> ^ ^iAtiH-iV,)lh _ 2iAt{H- 2Vo) 

h 
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Multiplying the above by Q' ^ and expanding the exponential terms in the Taylor series 
and finally rearranging the equation, one obtains: 


Q'Q ' t = 

h 


(2.76) 


Hence, 


< <S{t + At)\'^{t + At) > = 


»2VoAt/ft 


4( 


AtVo. 
h ‘ 




(2.77) 


It is clear that the error in the norm would grow exponentially. It is worth mentioning here 
that the use of a NIP shifts the eigenvalue spectrum to the complex plane and methods 
which are conditionally stable, and which depend on the eigenvalue spectrum become 
unstable under such conditions. 


2.4.3 The Split Operator (SO) Method 

In the application of the evolution operator, there is an error arising from the fact that 
the kinetic and potential energy operators do not commute. However, by splitting the 
time evolution operator symmetrically the commutator error is of third order only. Such 
an approach is known as the second order SO scheme [64]. It can be either potential 
referenced or kinetic referenced [76]. The potential referenced SO scheme is given by: 




(2.78) 


In the kinetic referenced SO scheme the exponetial containing the potential energy operator 
is symmetrically split and that containing the kinetic energy operator is sandwiched in 
between. The time evolution of the wave function in the potential referenced scheme is 
given by 


= Qm 


(2.79) 
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and 


Qq\ _ g-irAt/2ftg-iVAt/ftg-itAt/2R ^ g»rAt/2ftgiVAt/AgitAt/2?i 

= 1 (2.80) 


Hence the scheme is unitary, and the norm of the wave function is conserved. However, 
because of the non-commutability of the kinetic and the potential energy operators the 
scheme does not conserve energy. The scheme is unconditionally stable and it does not 
depend on the kinetic energy spectrum of the grid. However, the time step is selected based 
on the maximum value of the potential energy on the grid. Good results are obtained when 


At < 7r/(3AI4««) 


(2.81) 


where A14uu;(= ” ^tn) is the maximum excursion of the potential. This step size is 

chosen as to accommodate the entire energy level spectrtim of the system under investiga- 
tion [65]. 

In the presence of a NIP, Q' becomes: 


g/ _ g-irAt/2Ag-tVAt/ftg-VoAt/ftg-ttAt/2A 

The unitarity relation now becomes 

Q'Q' ^ = ^TAt/2h^-2VoAt/h^-itAt/2h 


(2.82) 


(2.83) 


In the absence of commutation errors the Icist two terms can be interchanged and the norm 
of the wave fimction at time t -|- At becomes: 

{'^{t + At)\^{t + At)) = ('$^(t)|Q'tQ|^(t)) 

= e-2^''^‘/''($(i)|^(t)) (2.84) 

It is clear that the norm decreases exponentially with increase in At and energy is no 


longer conserved. 
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2.4.4 The Chebyshev Polynomial (CP) Scheme 

Chebyshev polynomials are found to be superior to many other polynomials and are op- 
timal for a scalar function F{x) bounded in the interval [—1, 1]. So, a scalar function like 
e“® can be expressed in terms of these polynomials in the interval -1 < x < 1 as 

e- = f;(2 - <5„o)/n(a)i;(x) (2.85) 

n=0 

where 5„o is the Kroenecker delta and a = AEAt/2h. J„(q) are the modified Bessel 
functions of order n. Tn{x) are the Chebyshev polynomials of order n, calculated usi ng 
the recxirsion relation [233] 


(x) = 2xT„(x) - Tn_i (x) (2.86) 

with To(x) = 1 and T’i(x) = x. 

The evolution operator is a function of an operator. It can be shown (see Appendix 1) 
that a function of an operator can be expressed as a function of a scalar in the complete 
basis of the operator. So, the function of the operator can be approximated in the Cheby- 
shev series, provided the domain of the operator is confined to the interval [—1, 1] in which 
the Chebyshev polynomials are optimal. In case of a Hamiltonian which is self-adjoint, 
the eigenvalues lie on a real axis, and they can be positioned from -1 to 1 by renormaUzing 
the Hamiltonian as follows: 

^nonn = (2.87) 

where H = {Emax + E^in)/‘^ and AE = {E^ - Emin)- ^max and Emin liave been defined 
earher. In terms of this renormafized Hamiltonian, Hnorrm the evolution operator can be 
written as: 

g iHCii/h _ ioAnorm ^2 §§) 


The first term in the above expression is the phase shift due to the shift of the energy 



31 


scale. The second term is approximated by the Cheby^hev series as [69, 191, 196] 

OO 

^(2 - 5no) (2.89) 

n=0 

where $n(— i-Hnorm) are the complex Chebyshev polynomials of order n satisfying the re- 
cursion relation: 

^n+l ~ 2i^?norTn^n ^n— 1 (2.90) 

where $o = 1 and $i = —iHnonn- Therefore, the evolution of ^(t) in this scheme on a 
discrete grid is given by: 

+ At) = f;(2 - 5no) Jn(a)$n(-i#n.^)^(t) (2.91) 

n=0 

The number of terms to be used in the above expansion is estimated from the time - 
energy phase space volume a. In practice the number of terms used is slightly larger than 
this estimate for a good convergence. Since the evolution operator is expanded in a series 
of polynomials in the Chebyshev method by definition the scheme is not unitary. The 
deviation from the unitarity corresponds to the remainder term in the expansion. This 
deviation is used as an accuracy check of the scheme. The errors cire unif o rml y distributed 
in the bounded interval [191, 196]. Since Bessel functions show exponential convergence 
for n > a, the error is usually very small. 

The instability of the CP scheme in presence of a NIP was first noticed by Mowrey 
[234, 185]. It has been realised subsequently [235, 196, 236, 237] that the instability of the 
scheme is caused by the NIP. The main point to note here is that the Hamiltonian ceases 
to be Hermitian once the NIP is introduced. The eigenvalue spectnim of the Hamiltonian 
shifts to the complex energy plane and the renormalization of the Hamiltonian introduced 
above becomes invalid. 

Kosloff and coworkers [238, 196] have proposed the Newton interpolating polynomials 
in which the interpolation points are predetermined on a boimdsuy curve which encloses 
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the spectrum of the Hamiltonian. This is done by conformal mapping. For a Hermitian 
Hamiltonian the interpolation points become the zeroes of the Chebyshev polynomial, 
which is not the case for a non-Hermitian Hamiltonian. Huang et al [235] proposed the 
use of generalized Faber polynomials of which Chebyshev polynomials constitute a special 
case. The latter are generated by a one-to-one conformal mapping of the exterior of a disk 
to the exterior of a simple closed curve, L^. In case of the Chebyshev polynomials, 
becomes an ellipse. The renormalization of the Hamiltonian is done in a way as to account 
for the shift of its eigenvalue spectrum to the complex plane in the non-Hermitian case. 
These authors used the renormalization 

- ^H-iS , , 

-Hnorm = 2 ■ (2.92) 

where 27 > AE. The eigenvalues of Hnom are no longer bounded by one but axe mapped 
conformally from a unit disk. 27 equals the sum of the major and the minor axis. H in 
Eq. (2.92) is given by 

H = 1/2 [Em^ + Emin] “ (2.93) 

The evolution operator becomes 

^ ^-im/h £(_i)n(2 _ 6^)Jn{AEtJh)Tn{H^o,^) (2.94) 

n=0 

. Mandelstham et al [236] have made a simple analytic continuation of the Chebyshev poly- 
nomials while keeping their properties unaltered. Their method preserves the renormal- 
ization step as in the case of the Hermitian Hamiltonian in presence of the NIP. They have 

derived an expression for the time evolution operator as: 

N 

= S «n(*)Qn(-ffnorm; T) (2.95) 

n=0 

where an(t) = [(2 - 6r,o)e-^^^^\-l)^JniAEt/h)] and Qn{Hnorm , 7) are the analytic contin- 
uation of the Chebyshev polynomials satisfying the following recursion relation: 


C ^Qn— l(-ffnorm) T) "t" C^Qn+lC-HnormiT) — 2.Hn.ormQn(-HnorT7ij ‘T) 


(2.96) 
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with Qo(-Hnorm;7) = I and Qi{Hnorm'-,l) = e~^Hnorm- The Operator 7 is dimensionless and 
it defines the damping factor. This is connected to the NIP [P] in the expansion of the 
Green’s function by: 

r = Ai?[cos<^(l — cosh7) — isin<;6sinh7] (2.97) 

where <j) [= cos“^ phase. 7 is set to zero in the strong interaction region 

and it rises slowly as the asymptotic regions are approached. In this prescription, the 
NIPs are used externally as a damping function. Therefore, their complex nature does not 
interfere with the eigenvalue spectrum of the Hamiltonian. Another important point to 
be mentioned here is that the CP expansion method is conditionally stable because the 
scheme will be unstable if the energy range of the Hamiltonian is underestimated. The 
scheme does not conserve norm and energy and the resulting deviation is used as measure 
of the error [191, 196]. 

Very recently Neuhauser [237] has investigated the applicability of the NIPs in TIQM 
scattering. He proposes an anomaly-free very short range imaginary potentials which 
cover only 1-2 grid points. His method involves the computation of the full 5— matrix by 
diagonalization/inversion of the complex Hamiltonian. 

2.4.5 The Short Iterative Lanczos (SIL) Method 

The short iterative Lanczos (SIL) scheme, originally developed by Park and Light [72] 
computes the action of time evolution operator by forming a reduced subspace (Krylov 
space) of the Hamiltonian matrix. The orthonormal basis set, q^{j = 0, ...., N— 1) (Krylov 
basis set) which is spatially and temporally tailored, is generated from the initial vector 
9o = V’(O) follows [72, 239, 240]: 

Hqo — no^o + 0oQi (2.98) 


Hqj = ^j-\qj-i + cxjqj + /3jqj+i, j>l 


(2.99) 
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with 




( 2 . 100 ) 


and 


Pj-i = {Qj-i\H\qj) 

The Hamiltonian matrix becomes tridiagonal [71] in this reduced subspace: 

f OCq 00 0 0 ^ 

00 OL\ 0\ 0 * • • 0 

^ 0 01 a2 02 • • • 0 


( 2 - 101 ) 


( 2 . 102 ) 


0 0 0 • - • OCff^-2 0Nt-2 

\ 0 0 0 • • • 0Ni,-2 J 

It is of size Nl x Nl compared to N x N (where N is the total number of grid points, of 
the order of several hundred thousands for multidimensional scattering problems) of the 
original Hamiltonian. The Hermiticity of H is retained in is diagonalized by the 

unitary transformation matrix Z, the column of which contains its eigenvectors: 




(2.103) 


where is the diagonal matrix of eigenvalues and is the conjugate transpose of Z. 
The evolution operator can now be written as [75] 

U{At) = 

= (2.104) 


The time-evolved wave function then becomes 

= £'0-(At)4„ 

71=0 


(2.105) 
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The sununation here arises because the first beisis function, in the Krylov spacs, go is 
identical to ^(t). The above procedure is repeated iteratively N times. For small At (« 
0.1/s) typically 5-10 basis functions are sufficient to construct the tridiagonal matrix. This 
method is particularly useful for small At and since the scheme involves the exponential 
operation in the reduced subspace, it conserves norm as in the case of the SO scheme. It 
conserves energy also as the exponential operator commutes with the Hamiltonian unlike 
in the case of the SO scheme where it involves a splitting into kinetic and potentieil energy 
parts. The scheme is imconditionally stable. 

In the presence of a NIP the SH scheme is stable. But the basis vectors become 
nonorthogonal and one needs to resort to Gram-Schmidt orthogonalization at regular in- 
tervals. 

For large grid computations the method has additional disadvantage in terms of storage 
requirement. This is because the recursion vectors of size {Nz x N) need to be stored to 
evaluate the final time-evolved wave function. This problem can be circumvented if the 
recursion vectors are computed twice, once for constructing the tridiagonal matrix and then 
for use in the final equation to obtain the time-evolved wave function. In that case, the 
method becomes CPU intensive. With an increase in the mnnber of iterations the recursion 
vectors loose the orthogonality readily because of the truncation error introduced at each 
time stepv Error, starts btulding up when the recursion vectors span outside the reduced 
space. Park and Light [72] estimated this error by multiplying together the off-diagonal 
elements of the tridiagonal matrix: 

n A (2-106) 

where |gjVj.(Af)| is the magnitude of the first vector lying outside the Krylov space. 

A comparison of the various time evolution schemes in terms of their stability has been 
made by Leforestier et al [75]. We summarise our findings, with special emphasis on the 
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influence of the NIP on the various schemes. 

The SOD scheme is conditionally stable. It requires the Hamiltonian to be strictly 
Hermitian. It conserves norm of the wave function and energy. The error accumulates 
in the phase of the wave function. The scheme is unstable in presence of a NIP. The SO 
scheme is unconditionally stable. It conserves norm but the non-commutability of T and 
V” leads to non-conservation of energy. The exponential operation gives rise to a damping 
factor in presence of a NIP and it can be used successfully for damping of WP near the 
grid edges. Since the method does not depend upon the spectral range of the H amil tonian 
it is stable in presence of a NIP. The CP scheme is conditionally stable, and it does not 
conserve either norm or energy. The deviations axe. used eis a check on the accuracy of 
the method. The stability of the scheme depends strongly on the energy range of the 
Hamiltonian. The scheme becomes tmstable in presence of a NIP. This can be made stable 
by accounting for the shift of the eigenvalue spectrum of the H amil tonian to the complex 
plane in presence of the NIP. The SIL scheme also involves an exponential operation and is 
also unconditionally stable like the SO scheme. It conserves norm as well cis energy. This 
also results in the same damping factor eis the SO, in presence of a NIP. The accuracy 
of the scheme depends on the length of the time step. The additional complication of 
the non-orthogonality of the basis vectors in presence of a NIP can be taken care of by 
using the Gram-Schmidt orthogonalization scheme. The last point to be mentioned here 
is that the Hamiltonian becoming complex does not preserve the time-reversal symmetry 
(see Appendix 1) of the TDSE. This raises a fundamental question on the use of NIPs in 
solving the TDSE. We Avill show some numerical results to highlight this aspect in chapter 
3. 
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2.5 The Fourier Grid Hamiltonian Method 


The Fourier Grid Hamiltonian (FGH) method was proposed by Marston and Balint-Kurti 
[241] for computing bound state eigenvalues and eigenvectors using the variational princi- 
ple. In this approach the Hamiltonian matrix is discretized on a set of uniformly sampled 
points on the discrete Hilbert space and the forward and the reverse Fourier transforms are 
reduced to a sum over cosine functions. The diagonalization of the kinetic energy operator 
in the FGH method, as can be seen below, is a special case of the DVR method. 

For a one dimensional system, we can write 


/*oo 

= / {x\k) — {k\x')dk + V{x)S{x - x') 

J—oo 2771 

= r° + V(x)6(x - x') 

2x7-00 2m ^ ^ ^ ' 


The last equality is obtained by using 


(2.107) 


In finite space, the above equation can be written as 


(2.108) 


(Xi|.^Xj) 


1 A r^'“(ZAA:)2) _ V{xi)6i, 


±_ y- giiA;fe(x.-x,) ) >•' 

2x , 1 2771 




AA: + 


Ax 


= ^ + v(x,)4, 


(2.109) 


where Tj = ^ 2 ^^* Combining the negative and positive values of I one obtains 


iriy= ;^{|^E“s(i2>r(i-j)/JV)T, + V(xi)i„| (2.110) 

since To = 0 by definition. Any arbitrary state vector \<j> > can be written as a linear 
combination of the eigenvectors |xj > for the discrete grid points: 


1^ > = ^ |aJi > Ax < Xi\(f> > 

t 
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= ^ !®i > Aa;<?i(a:i) (2.111) 

i 

where 4>{xi) are the values of the wave function at the grid points and are the unknown 
coefficients to be evaluated by the Arariational method: 

p = < ^ 1 -^^ > 

<(^!^> 

_ Eij <f>(xiyAxHijAx6{xj) 

Ax'£i\<l>(xi)\^ 

E, ^ ' 

where = HijAx. Minimization of this energy with respect to the variational parameters 
<l){xi) yields the set of secular equations: 

•£ K - “ (2,113) 

j 

The nontrivial solution for the <f>{xj)s are obtained if 

- ESij) = 0 ( 2 . 114 ) 

This is the characteristic equation having A roots for E. 

In the numerical calculation we have stented by defining an eigenvector matrix for 
(j) containing A columns. The matrix is composed entirely of zeroes except for a single 
element of unity which occurs at the A*^ row. The matrix is computed using these 
eigenvector matrices in which a forward and a reverse FFT are performed for computing 
the action of the kinetic energy operator on <f). This Hamiltonian matrix is then 
diagonalized using the standeird matrix diagonalization routine. The diagonal elements 
of the resulting Hamiltonian matrix are the A eigenvalues and the columns of (ft are the 
corresponding orthonormal eigenvectors of The resulting eigenvectors are renormalised 
using a standard integration procedure. 
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2.6 The Spectral Quantization Method 


The spectral quantization method (SQM) was proposed by Feit et al [65, 242] for computing 
the bound state eigenvalues and eigenfunctions by solving the TDSE. In this method a 
WP (usually a stationary Gaussian wave packet (GWP)), ^(0), is located initially in the 
interaction region of the PES. The WP at time f, is obtained by evolving $(0) in 

space (by the FFT method) and in time (by the SO method). The temporal autocorrelation 
function 

C{t) = (^(0)|^(*)) (2.115) 


is computed at various time intervals and Fourier transformed to generate the power spec- 
trum: 


2 

1(E) = jT C(t)e^l^dt 


(2.116) 


1(E) is known as the spectral intensity. The peaks in the power spectrum arise from the 
bound and/or quaisibound states of the system and the energy values corresponding to the 
peak maxima are their eigenvalues. The energy resolution of the spectrum is decided by: 

27r/i 




(2.117) 


where T(= Nt x At) is the total time up to which the WP is evolved. Fomier transform 
of C(t) leads to a superposition of 5-functions [65]: 


1(E) = 2 - ^n) (2.118) 

where = Wn, are the relative weights of the states. This 5-function behavior 

implies an i nfini te record of the autocorrelation function data. In order to correct for the 
finite time computations, the autocorrelation function is multiplied at each time step by 
the normalized Hanning window fimction w(t)/T where 

w(t) = 1 — cos 0 ^ ^ ^ 

if t > T 


0 


(2.119) 
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before Fourier transformation. The use of the window function improves the resolution of 
the spectral peaks. In the absence of the window function, the spectral peaks resemble the 
function [65], ^nd their soulders overlap with the nearby peaks. The accuracy 

of the eigenvalues is improved by fitting I{E) to a Hneshape function L\{E—En) assuming 
a single-line fit: 

I{E) = WME - (2-120) 


where L\{E — JE?„) is given by [65]: 


Li{E-Er,) = exp[i{E - En)t]w{t)dt 

exp[i{E - En)I] - 1 


_1 

2 

1 

■^2 


exp{[i{E — En)T-h 2tt]} - 1 


i[{E-E^)T+2'K] 
exp{[i(J5' — En)T— 2 tt\} — 1 


( 2 . 121 ) 


i[(E-F;n)T-27r] 

T his single line fit is described in detail by Feit et al [243]. 

Once the eigenvalues are computed, the eigenfunctions can be obtained by projecting 
the time evolved wave packet at the desired eigenstate En~ 

rT 


~ 77, / ^{t)w{t)exp{iEnt)dt 
1 Jo 


( 2 . 122 ) 


The eigenfimctions obtained in this way are renormalised by using either the Simpson’s 
rule or by some other quadrature scheme. It is worth mentioning at this point that the 
eigenfunctions obtained by this method may have an error in their phase, arising from 
the use of the second order SO scheme. Our experience shows that these eigenfunctions 
lead to inaccuracies in the reaction probability calculation. The phase error vanishes 
when the probability density of the eigenfunction is computed and there is no problem in 
characterizing the states based on the nodal pattern. We have verified this by computing 
the boimd state eigenfunctions for Hj using an extended Rydberg potential, both by the 
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FGH [241] method and the SQM [65]. We find that the eigenfunctions obtained by the 
two methods do not differ in terms of their probability density patterns but do show 
the difference when they are used in reactive scattering calculations. The eigenfunctions 
obtained by the FGH method are found to be quantitatively accurate. Bandrauk and Shen 
[244] proposed the exponential SO method which is accurate up to third order: 

g-i(t+V 0 t/fi ^ g 7 Ar/ 2 g 7 AVg(l- 7 )Af/ 2 g(l- 2 T)AVg(l- 7 )Af/ 2 g 7 AVg 7 AT '/2 (2.123) 

where 7 = ^ 2 - 21 / 3 ) ^ ~ Their analysis shows that use of this exponential SO in 

place of the second order SO removes the phase error in the eigenfunctions. The SQM 
is used to compute the resonances in the present thesis and the relevant computational 
details eue given in chapter 3. 

2.7 Scattering Analysis 

2.7.1 Time-Energy Mapping of the Wave Packet Flux 

In scattering studies the reactants are initially represented in terms of an appropriate 
wave packet located far out in the asymptotic reagent channel where the the interaction 
potential is almost zero. The configuration space is divided into reagent, interaction and 
product regions and the vibrational state-selected energy resolved reaction probabilities 
( Ff{E!) ) are calculated by computing the energy resolved flux of the WP across a divid- 
ing line in the product channel. This strategy was developed initially by Neuhauser and 
Baer [155, 156, 164] and later on was adapted by others [166, 167, 245, 246]. In collinear 
(A,BC) type collisions the initial WP ($(i2,r,f = 0) in mass-scaled Jacobi coordinates R 
(A, BC translation) and r (BC vibration) consists of a product of a minimum uncertainty 
Gaussian wave packet (GWP), F{R), corresponding to (A,BC) translation and the vibra- 
tional eigenfunction, ^v(r), corresponding to the vibrational state v of the BC molecule; 


^R,r,t = 0) = F{R)Ur) 


(2,124) 
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where 


FXR) = (2'irS^y^^exp 


{R-Roy 


— ikoR 


(2.125) 


Here 6 is the width parameter, Rq and ko denote the location of the TnaviiniiTn in the 
coordinate and the momentum space respectively. The vibrational eigenfunctions, 
of BC are nmnerically computed using the FGH method [241]. The reaction probability, 
lf(£^ , is computed from 



U(R,r,,E) 


mR,r,,E) 

dr 


(2.126) 


where p = [m^mBmc/(m^ + tub + is the three body reduced mass and the quan- 

tity within the angular bracket is the energy-resolved flirx of the WP (computed along the 
dividing line at r = r/ in the product channel), which is integrated over the entire range of 
R. ^(i^, Tj, E!) is the energy-resolved wave function normalized with respect to the initial 
translational energy distribution of the WP, i.e.. 


^(R,rj,E)=^(R,rj,E)/aE (2.127) 

$(i?,r/, E) is obtedned by Fourier transforming the time evolved WP, $(JJ,r,t), edong the 
dividing line at r = r/ : 

^(R,ri,E) = —^ [ ^(R,r,t)exp(iEt/h)dt\r=rf (2.128) 

vStt J—oo 

as in Eq. (2.127) is the weight of the energy component, E, contained in the initial 
translational WP and is evaluated by 

as = (2.129) 

where 


F{k) = (25^/7r)^/^exp [-^^(fc - ko)^ + i{k - fco)i2o] 


(2.131) 
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The factor {n/hkY^'^ in Eq.(2.130) accounts for the energy normalization of ^{R,ri,E). 
Incorporating the result of Eqs. (2.127) - (2.131) in Eq. (2.126), we obteiin 



^{R,ri,E) 


d^{R,ri,E) 


dr 


(2.132) 


with k = \j2ix{E— Ey)/h, and E^ the initial vibrational energy of the BC molecule. 

In actual computation, ^{R, rj, E) and its derivative with respect to r are computed at 
each time step along r = rjhy evaluating the continuous Fourier transform (Eq. (2.128)) 
as a discrete sum as follows: 


At °° 

#(i?, rj, E) = $(i?, r, n A t)e 

v27r„=o 


iEnAt/hi 


r/, E)| ^ At d9{R,r,n At) iEjuxt/hx 

EiZ ^ ^ l»^ 

v27r„^ 


(2.133) 


(2.134) 


Such an approach reduces the storage requirement of the computer to the size of the R grid 
times the number of energy values at which Ff{E) is sought. Once the contributions from 
each time step are added coherently to the values from the previous step in Eq. (2.133) 
and (2.134), the WP is damped out with the help of either a NIP or a masking function. 


2.7.2 Time-Energy Projection 


In this method the energy resolved state-to-state reaction probabilities (Pv,,y(E)) are com- 
puted from the energy resolved matrix elements for the transition from a specified 
initial reactant level (u) to a specific product level {v'). The 5^,,/ matrix elements in this 
method are computed by solving the TDSE. The method was proposed by Balint-Kurti 
and coworkers [162, 227]. These authors have applied the method initially for F -h H 2 
[162] and later for Li + FH [247, 231] reactive scattering. For the collinear (A,BC) type 
of collisions the '9{R,r,t = 0) is prepeired in the same way as mentioned in the previous 
subsection. This WP is time evolved and the resulting '^{R,r,t) is converted to ^ in 
product coordinates (i2',r'). One may use a Fourier type of interpolation scheme [164] in 
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which ^(i2, r, t) is first transformed to the momentum domain '‘^^(pR,Pr, f) through an FFT 
route and then projected on to the product coordinates through an inverse FFT : 


^(iT, r', t) (2.135) 

PRJPr 


where 

*(PB,Pr,f) = (2.136) 

Nr and Nr are the nmnber of grid points along R and r respectively. This interpolated 
WP at a fixed value of the product scattering coordinate i2j, ^{Rfi, r', t), is projected onto 
the specific vibrational state (t/) of the product molecule to extract the contribution of 
that particular state in This yields a time-dependent coefficient at each time: 


max 

CAt)= I <f>Ar'MKr',t)dr 


(2.137) 


where are generated in the same way as (^vi^) by the FGH method [241]. The 

Fourier transform of (f ) results in energy dependent amplitude 

A^(E) = ^ r e'^^^CAt)dt (2.138) 

ZTT Jo 

from which the P^^^{E) values are calculated: 


P,y(£0 = |5'.y(P)P = -Kk^ 




F{K) 


(2.139) 


In this method the interpolation of the WP takes much of the CPU time but, storage 
is no longer a problem. The CPU requirement can be reduced if the WP is interpolated 
only once when it reaches the interaction region and then the dynamics is monitored in the 
product grid. This also reduces the error arising from interpolation, as the entire process 
requires only one interpolation and that too in a compact form of the WP. 


2.7.3 The M0ller Wave Operator Formalism 

Viswanathan et al [139] have described computation of the 5-matrix elements through 
a time-dependent wave packet (TDWP) calculation employing M0ller wave operators. 
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Weeks cind Taimor have [140, 141, 142] adapted their technique to compute the scattering 
attributes. The Mpller operators are defined as: 

ai = (2.140) 

where 7 refers to the channel index. ^{= “IlR) + h'’'(r)) and H{= + V) are the 

asymptotic and the full interaction Hamiltonians, hP'{r) is the diatomic Hamiltonian in 
the channel and V' is the interaction potential. These Mpller operators are used to 
prepare the appropriate scattering states of reactants (products). For the reactant packet 
one initially starts with ^(i2, r, f = 0) centered at the interaction region of the PES {t — 0) 
and operates with 0" on it to prepare the appropriate reactant channel (a) WP, 

The action of on ^ ( JZ, r, f = 0) is to propagate it backward in time under the action 
of the asymptotic Hamiltonian (Hq) and then to propagate it forward in time again to 
the interaction region imder the action of the full interaction Hamiltonian (H). In actual 
computation, the reactant packet is back propagated for a length of time t = — ri which 
takes it far out in the asymptotic reagent channel and then propagated forweird in time 
to the interaction region. This forw2U-d propagation is stopped before the WP starts 
bifurcating. This whole process is called the reactant channel dynamics. A similar exercise 
is undertaken to prepare the appropriate product channel (/3) WP, 'St, and this is called 
the product channel dynamics. In this case the initial WP is written in terms of the 
product Jacobi coordinates, ^'(i2',r',f = 0) and located in the interaction region of the 
PES and Qt is acted on it to obtain The initial propagation under the action of the 
appropriate asymptotic product channel Hamiltonian {H^) is carried out forward in time 
(f = Ti) on the product grid and brought back to the interaction region by propagating 
backward in time and the propagation is stopped before the WP starts bifurcating. Once 
these M0ller states are obtained the dynamics is followed on a common grid. 

The common grid can be obtained either by interpolating the reactant chmnel packet 
to the product grid and product channel packet to the reactant grid or by interpolating 
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both the channel packets to a common grid. A Fourier type of interpolation [164] as 
mentioned earlier can be used for this purpose. The reactant packet is then propagated 
forward in time from t = 0 to i = r and its correlation function with the product channel 
packet (frozen at t=0) is computed: 

(2.141) 

Similarly, the product channel packet is propagated backward in time from t = 0 to t = — t 
and its correlation function with the reactant channel packet (frozen at t=0) is computed: 

Cap{t) = (2.142) 


When these two parts are combined together it results in the correlation function between 
the reactant and product channel packets over the entire time range from — r to t. The 
symmetric partitioning is possible because of the symmetric form of the correlation function 
= Ca 0 [140], The time r corresponds to the value of time when the correlation 
function becomes zero. The components of the WP which reach the grid edges during the 
correlation function computation are absorbed using a NIP. The correlation function is 
Fourier transformed to compute the energy resolved S'-matrix elements: 




(2.143) 


where the F terms are computed using the Eq (2.125) and fcX = . The 

state-to-state reaction probabilities are computed from the S'-matrix elements as: 


P^(F0 = (2.144) 

An alternative way is to start initially with a rectangular bond coordinate grid (r^B> ^bc)- 
The reactant and product cannel WPs are then prepared by switching over to the respective 
grid through interpolation and then the channel WPs are transformed back to the bond 
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coordinate grid where the rest of the dynamics is followed in the same way as mentioned 
above with the kinetically coupled Hamiltonian [248]: 


^ab dr%c 2/iB drAB ^Bc 


(2.145) 


This approach exploits the collinear constraint on the dynamics and is better because 
the grid spacings remain uniform. Jackie and Meyer [143] have recently applied this 
correlation function approach of computing the ^-matrix elements in the bond coordinates 
to the collinear H + H 2 reactive collision in the frame work of multi-configuration time- 
dependent Hartree approach and obtained results which are in excellent agreement with the 
time-independent results. We have developed a code to compute the state-to-state reaction 
probabilities for the collinear He + Hj collisions using the M0ller operator formalism and 
we are in the process of debugging it. 



Chapter 3 


Transition State Resonances and 
Collinear He, (HD+/DH+) 
Dynamics 

3*1 Introduction 

Resonances in (He,H^) collisions have received considerable attention from several resemch 
groups during the past two decades. Quantxim mechanical studies of the collinear reaction 

He + Hi-(t-) HeH+(t;') + H (Rl) 

by Kouri and Baer [249] revealed several oscillations in the reaction probability, P^, versus 
energy (E) plots in the range 0.93-1.4 eV on a diatomics-in-molecules (DIM) potential- 
energy surface [250]. Subsequent calculations by Adams [251] using the integral equation 
approach on the same PES, on a finer mesh of E, in addition to reproducing the results 
of Kouri and Baer, revealed additional resonances. Chapman and Hayes [252] discovered 
several oscillations in the vibrationally inelastic transition probability versus energy plot 
in the energy range 0.4-0.65 eV, well below the reaction threshold. They characterized 
one of the oscillations to be a shape (open) resonance and the others as of Feshbach 
(compound) type. Joseph and Sathyamurthy [253] carried out R-matrix calculations for 
the system in collinear geometry on an accurate ab initio PES due to McLaughlin and 
Thompson [254]. Their results also revealed several oscillations in the P^(E) curves and 
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most of them were interpreted in terms of quasi-bound states supported by vibrational 
adiabatic potentials [255]. Kress et al [256] adopted the adiabatically adjusted principal 
axis, hyperspherical formulation of Pack et al [257] for quantum reactive scattering in three 
dimensions (3D) eind reported the results jfor total angular momentum, J = 0, collision 
on the McLaughlin-Thompson-Joseph-Sathyamurthy [258] (MTJS) PES. Their results re- 
vealed ntimerous resonances with life-times of ~ 1 ps or more. Similar calculations were 
carried out by Lepetit and Launay [259] on the same PES and they reported some of the 
resonances having life-times ~ 3 ps. Results of Zhang et al [260] for 3D, J = 0 collisions 
on the DIM PES also showed several resonances. 

Stroud et al [261] had carried out a time-dependent quantum mechanical (TDQM) 
investigation of the collinear reaction on a spline-fitted ab initio PES. Balakrishnan and 
Sathyamurthy [262] carried out a TDQM calculation for the collinear collision on the 
MTJS PES and they also obtained several resonances. These were related to chaotic 
dynamics in classical calculations [263, 264] and turbulence in the flux patterns obtained 
in the quantal calculation. They also reported [265] the preliminary results of a 3D, 
J = 0 quantal investigation of the reaction (Rl). Although they did not report energy 
resolved values, the probability density plots obtained by them were highly structured, 
suggesting the existence of dynamical resonances for the system. Mandelshtam et al [266] 
have, traced the connection between the quantum mechanical resonances and the classical' 
chaotic scattering in three dimensions. Marston [267] have interpreted the resonances in 
the chaotic region in terms of periodic orbits of the collinear HeH^ system. More recently, 
we have extracted the signature of quantum chaos in collinear (He, H 2 ) collisions [268]. 

Recently Sakimoto and Onda [269] obtained the resonances in the system for the 
collinear geometry on MTJS and DIM PESs using a time independent quantal calcula- 
tion. While the positions of resonances computed by them were in agreement with those 
of Sathyamurthy et al [255] the magnitudes of P^ were at variance. More recent calcula- 
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tions through [167] time-energy mapping of the reactive flux of the WP )rield I^{E) values 
which are in quantitative agreement with those of Sakimoto and Onda. 

In the present chapter, we analyse the dynamical resonances in collinear (He,Hj) col- 
lisions and its isotopic variants when either one of the two H atoms is replaced by D, in 
terms of their transition state spectrum. Recently, Sadeghi and Skodje [97, 98] showed 
that they could obtain thirteen high energy resonances for collinear (H,H 2 ) system by a 
correlation function approach. They solved the TDSE for the system by the spectral quan- 
tization method [64, 65, 242] to generate the transition state spectrum for the dynamics 
on the double-many-body-expansion (DMBE) potentiail-energy surface. We have extended 
their approach to a study of collinear (He,!!?) dynamics and its isotopic variants. We 
have also calculated the vibrational (u) state-selected energy resolved reaction probabili- 
ties ( (E)) for collinear (He, HD''') and (He, DH'*') collisions by the time-energy mapping 
of the reactive flux [155, 156, 164, 167] of the WP as described in Sec.2.7.1. 

3.2 Computational Details 

3.2.1 Computation of the Transition Spectrum 


We have used the SQM outlined in sec. 2.6 for the computation of the transition state 
spectrum. The interaction Hamiltonian for the collinear A + BC system in mass-scaled 
Jacobi coordinates {R,r) is given by 




+ V(B,r) 


(3.1) 


2n [dB!^ dr^j 

where // is the three-particle reduced mass. The MTJS potential-energy surface, used for 
Y(R,r), has a well of depth 0.32 eV in the interaction region. The TDSE is solved on a 
256x256 grid in (R,r) space with the origin at (2.2,0.4) a.u. and AR = Ar = 0.05 a.u. 

The initial wave function, ^(i?,r,t = 0) was taken as an even parity Gaussian wave 
packet (GWP) in terms of A-B and B-C bond coordinates: 



^(iJ, r,t = 0) = N exp 


r i^AB - r%)^ 


262 


(f'BC ~ fBc)^ 
262 


(3.2) 


where N is the normalisation constant and ( r^, ) specifies the initial location of the 

WP. R and r are related to and rgc through: 


''"AB 


f'BC 


dR- 


rBC 

2 


d~^ r 


(3.3) 


where d is the scaling parameter, given by the square root of the ratio of the reduced mass 
of BC to the three body reduced mass p and 6 is the width parameter. The average energy 
of the WP, dependent on its location and its width, is obtained from 


(^(fl,r,t = 0)W$(fl.r,t = 0)) 

' ' {4'(iJ,r,t = 0)|4'(iJ,r,e = 0)) 

The value of Af used in the SO method for time evolution was selected on the basis of 
the desired width of the energy spectrum. The maximum width in energy that could be 
obtained is given by 


A Emax — 


Tvh 

Af 


(3.5) 


We have carried out six sets of calculations for the HeHj and three sets each for the HeHD'*' 
and HeDH"^ by varying the position of the initial WP in the interaction region. We have 
used At = 0.161625 fe for each set except for the sixth set of the first system, in which Af 
was taken to be 0.10775 £s. The lengths of the time step 2 ire selected in accordance with 
the Eq. (2.81). The grid parameters for the different calculations are given in Table 3.1. 

We have used the NIP of Neuhauser et al [200, 201] in both reactant and product 
ch ann els in the asymptotic region to suppress the reflection of the WP from the grid edges 


— ' — irwTlil 
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at a later time. The form of the potential tised in the calculation is the linear ramp: 


Vi{X) = 



X-Xu 

Xii-Xu 


i£Xu < X < X 21 
otherwise 


(3.6) 


where Vj(X) is the absorbing potential applied at the channel coordinate X, Vq is a real 
constEint specifying the height of the potential, Xi/ is the point where the absorbing bound- 
ary starts and X 2 / is the end of the boundary. X 21 — Xu = AXj is thus the width of the 
boundary. The height and the width of the absorbing boimdary were chosen to satisfy the 
inequality condition [200, 201, 270, 271]; 


< 


AXi^M^E^ 


The left-hand inequality in the above equation guarantees complete absorption of the flux 
and the right-hand inequality guarantees that no reflection takes place due to interaction of 
the wave function with the imaginary potential. Choice of Vq = 0.25 eV and AXi = 2.75 
a.u. was found to give 99% absorption. We have also used a quadratic NIP used by 
Seideman and Miller [137]: 




ifXu < X < Xn 
otherwise 


In this case values of Iq = 0.25 eV and AXj = 2.25 a.u. were found to damp the WP 
more smoothly than the linear potential. 

The autocorrelation function for each initial choice of the WP was computed at each 
time step, multiplied by the normalized Hanning window function and Fourier transformed 
to the energy dommn to generate the transition state spectrum as described in sec. 2.6. 

We have time evolved for a total of 32768 time-steps in each set of our calculations. 
That corresponds to a total time of 5.29 ps. This gives an energy resolution of the order of 
7.8 X 10“'^ eV (Eq. (2.117)). This duration of time evolution was found to be suflicient to 
get most of the wave function out of the interaction zone, thus bringing the autocorrelation 
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No. 

r“43(a.u.) 


<E> (eV) 

HeH^ HeHD^ HeDH+ 

1 

2.133 

3.399 

0.6928 

0.6632 

0.6332 

2 

2.462 

3.567 

0.9924 

0.9628 

0.9332 

3 

2.533 

3.860 

1.2159 

- 

- 

4 

2.874 

4.529 

1.7693 

- 

- 

5 

3.345 

5.283 

2.2816 

2.2520 

2.2224 

6 

3.846 

6.414 

2.6834 

- 

- 


Table 3.1: Parameters for different choices of initial wave functions used in time evolution. 


function down neairly to zero. The transition state spectrum, 1(E), exhibits several sharp 
maxima, positions of which give the eigenvalues of the resonance states. Since the system 
possesses a large ntunber of closely spaced scattering resonances, some overlapping contri- 
bution remains even after multiplication by the window function. To make an accurate 
estimate of the eigenvalues, we have fitted 1(E) to the lineshape function (Eq. (2.121)). 

The eigenfunction at any given eigenenergy was computed through Eq. (2.122). The 
life-time of some of the well resolved states are computed by a stable nonlinear le 2 ist-squares 
fitting of the spectral peaks to a Lorentzian function: 




^max V 2 


Tn ^2 




{E-E^r+e-fY 


(3.9) 


where gives the width of the n‘^ resonance state, and the life-time for the n‘^ state is 
related to its width by 



(3,10) 


3.2.2 Computation of Reaction Probabilities 

The Py{E) values are computed by the time-energy mapping of the reactive flux of the 
WP (Sec. 2.7.1) for the collinear (He,HD'^) and (He,DH''‘) collisions. The vibrational 
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eigenfunctions (I>v{t) of HD'*' are numerically computed by the FGH method (Sec. 2.5) 
using the diatomic Hamiltonian: 


/i = 


2mhz)+ 


dr^ 


+ V(r) 


(3.11) 


where V(r) is the extended Rydberg potential used in the MTJS PES [258]: 


X^r) = — Z)e j^l + o,\P ■!“ U2P^ “b <^3P^ 


X e 


-aip 


p = r 


(3.12) 


Here is the equilibrium bond dissociation energy of Ht and ai, a 2 , are the parameters 
reported elsewhere [272]. We have used the FFT algorithm for the spatiaJ. evolution of the 
WP and the SO method for the temporal evolution. We have used a masking function 
activated outside the dividing line in the product channel and in the asymptotic reactant 
channel in order to damp out the WP from grid edges: 

(3.13) 

where Xmax is the maximum length of the grid along the coordinate as, Xf is the point 
at which the function hcis the maximum value (=1) and Ax is the width over which the 
function decays from 1 to 0. f{xi) values are multiplied with the final time evolved WP in 
each channel which results in a sin" masking. The WP is time evolved for a total of 773 
£s which is found to be sufficient to obtain converged Ff{E) values. The grid parameters 
used in this calculation are given in Table 3.2. 


/(®i) 


Sin 


^ (^mai ®t) 


Ax 


3.3 Results and Discussion 

3.3.1 Bound states and Resonances in HeHj 

An initial WP, (No.2) with {E)= 0.9924 eV, was located at = (2.462,3.568) 

a.u. and time evolved as described in section 3.2. Snapshots of the time-evolved wave 
function, at different time intervals are shown in Fig.3.1(a-1), in the form of probability 
density contours, superimposed on the potential-energy contours. In about 8 fe, the WP 
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Nr 

256 

Nr 

256 

(Rmin, iZmax, Ai?)(a.U.) 

(2.2, 14.95, 0.05) 


(0.4, 13.15, 0.05) 

Ro (a.u.) 

10.0 

5(a.u.) 

0.25 

■ At(fe) 

0.19 

rj (a.u.) 

8.4 

JfZnuuib ( a. U. ) 

12.95 

^ma»ife(^'^-) 

11.15 


Table 3.2: Grid parameters used in the scattering calctilation. 


spreads along the asymmetric stretch direction and in another 8 fe, it develops interesting 
structures. These structures are amplified with further time evolution and part of the 
wave fimction spreads in the (He,H^) channel and part in the (HeH'*',H) channel. Since 
we have used NIPs (eq. (3.6)) in both channels we do not see any back reflection of these 
outgoing waves firom the grid edges. As discussed by Balakrishnan and Sathyamurthy [262] 
substantial portion of the WP lingers on in the interaction region, revealing characteristic 
nodal structiires which indicate high vibrational excitation of Hj and also vibrational 
excitation of the HeHj complex. Even after 5.296 ps, there is 0.05% of the WP remaining 
in the interaction region indicating that some portion of the WP has a lifetime greater 
than 5.3 ps. 

We have computed the autocorrelation function (Eq. (2.115)) during the entire time 
evolution and it is clear firom Fig. 3.2 that a sharp decline in the first few fe of evolution 
is followed by a large number of recurrences over a long time (5 ps). The inset in Fig. 3.2 
shows clearly the magnitude of recurrences in the ps time scale. These recurrences can be 
interpreted in terms of a number of symmetric and asynametric stretch orbits for the 










R (au) R (au) 

Figure 3.1: Probability density contours of the initiafGaussian wave packet (No.2) with 
<E> = 0.9924 eV and width, S = 0.25 a.u. at different times t (indicated inside each box) 
showing its evolution on the (256x256) potential grid. The total probability density of the 
WP retained on the grid at different times are: 1, 1, 0.99, 0.99, 0.77, 0.035, 0.01, 6.5x10"^ 
and 5.3x10“"^ in a, b, c, d, e, f, g, k and 1 respectively. 
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system and the time interval between the recurrences can be connected to the periods of 
such classical orbits [273]. The natural oscillation of, |C(t)|, persists up to 4.85 ps. 

The spectral intensity, 1(E), (Eq. (2.116)) obtained by a half-Fourier transform of the 
above autocorrelation data, in the negative energy range is shown in Fig. 3.3. There are 
two sharp intense peaks corresponding to the two botmd states (-0.1281 eV and -0.0328 
eV). These values are in agreement with those computed by Sathyammthy et al [255], 

The spectrum in the energy range 0-3 eV is shown in Fig. 3.4 in the form of a logio/(E0 
plot. The intense peaks therein would correspond to quasi-bound states and hence to 
resonances. The spectrum shows a decaying pattern in the energy range beyond 2.4 eV, 
presumably due to increased noise arising from the fact that we have used a cut-off value 
of 3 eV for the potential on the grid. Such spectra are expected to be reliable up to 80% 
of the cut-off potential. The spectrum shows a complicated pattern due to the existence of 
a large number of closely spaced quasibound states for HeHt • Overlapping contributions 
to the spectral lines from the neighbouring states persist even after filtering, as indicated 
by the presence of shoulders on some of the spectral peaks. 

We have superimposed, in Fig. 3.5, the 1(E) values in the range, E=0.90-1.30 eV on the 
reaction probability {F^{E)) values for the ground vibrational state of Hi", reproduced 
from Ref. [255]. It is clear that some of the petiks in the P^(E) results must be actually 
made up of several narrower peaks. Alternatively, R-matrix calculations on a finer mesh 
of E-values would have revealed additional resonances. 

It is important to point out that the intensities of the different spectral peaks are 
dependent on the choice of location of the initial wave function. The latter decides the 
(E) and hence the contributions from different quasibound states. A comparison of the 
power spectrum obtained from the different sets of calculations reveals the commonality 
in the resulting spectra. It is also clear from these spectra that there is a sharp resonance 
near the threshold for each v and t/ state of H 2 and HeH"*" respectively. They are all 
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Figure 3.3: Transition state spectrum in the energy range -0.15-0.0 eV computed with the 
initial Gaussian WP No.2. Two distinct peaks in the spectrum correspond to the two true 
bound states of the HeH^ complex. Intensity, 1(E) is plotted in arbitrary units. 




log,oi(E) 



gure 3.4: Power spectrum in the energy range 0.0-3.0 eV computed with the initial 
aussian WP No.2. 
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listed along with the vibrational eigenvalues for Hj and HeH^, in Table 3.3. It is worth 
re-iterating that we also observe a number of resonances which do not correspond to the 
thresholds. 

Results obtained using the WP with (^?)=2.6834 eV is particularly interesting. As 
can be seen from Fig. 3.6, the decay of the autocorrelation function and the spectral 
pattern, shown as an inset are much simpler than in Fig. 3.2 and 3.3 respectively. The 
power spectrum stops abruptly aroimd 2.78 eV (the dissociation threshold, and decays 
smoothly thereafter. It is also clear from Fig. 3.6 that there is a large number of resonances, 
just below This is to be contrasted with the lack of resonances at higher energies 
reported by Sakimoto and Onda [269]. This could be partly due to the fact that the 
correlation function predicts the resonances but it does not give any indication of the 
magnitude of the oscillations in the P^(£?) curve. 

Since our calculations reveal a large number of resonances corresponding to the quasi- 
bound states of the HeH^ transient we have analysed some of the eigenfunctions. First 
we show the probability density contours for the two true bound states in Fig. 3.7(a) and 
3.7(b). Then we show the contour plots of for the two low-lying quasi-bound 

states in Fig. 3.7(c) and Fig. 3.7(d). The former is dominated by motion along R while 
the latter is dominated by motion along r. Such states can be characterized by a set of 
two quantum numbers: ^h+)- Thus the eigenfunctions depicted in Fig. 3.7(a-d) 

would correspond to (0,0), (1,0), (2,0) and (0,1) in that order. 

The wave functions plotted in Fig. 3.7(e-h) show an interesting behavior that is to be 
expected of those for quasibound states. With increase in R, the nodal structure along 
r decreases. For example, in Fig. 3.7(e), at small R, the wave function has one node 
in T direction but with increase in R the node along r ceases to exist. Similarly in Fig. 
3.7(f), the wave function exhibits one node along r at short R values and no node along 
r at large R. In Fig. 3.7(g) and 3.7(h) the number of nodes along r decreases from 2 to 
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En(eV) 

Et^CeV)* 

En(eV) 


0.1458 

0.1433 

0.9522 

0.9543 

0.4147 

0.4182 

1.3138 

1.3153 

0.6751 

0.6751 

1.6369 

1.6362 

0.9140 

0.9161 

1.9193 

1.9181 

1.1456 

1.1422 

2.1634 

2.1616 

1.3540 

1.3540 

2.3661 

2.3668 

1.5525 

1.5520 



1.7494 

1.7364 



1.9074 

1.9074 



2.0647 

2.0648 



2.1933 

2.2086 



2.3386 

2.3385 




Table 3.3: Eigenvalues (En) of the resonances near the threshold of various vibrational 
channels for He+H^ (v) — > HeH‘''(t/)+H reaction. ♦ Results of Ref. [255]. 


1 to 0 with increase in R. Therefore, it becomes difficult to assign clear-cut quantum 
numbers based on local modes for those eigenfunctions. With ftirther increase in energy, 
the motion can no longer be characterised in terms of R or r only. It acquires the charac- 
teristics of hyperspherical modes [274], as can be seen from Fig. 3.8(a) for the quasiboimd 
state at R„=0.9661 eV, slightly above the reaction threshold. An attempt to compute the 
resonant periodic orbit (RPO) [275] at that energy led to locating one at Erk)= 0-8861 
eV, the difference arising presumably due to lack of zero-point-energy in the classical orbit 
[276]. 

The resonance eigenfunction shown in Fig. 3.8(b) corresponds to the threshold of r/=l 
and u=4 and therefore shows clearly one node in the HeH'*' channel and 4 nodes in the 
channel indicating that it is indeed a predissodative state of =1 and v=4. 

As is expected of quantum states, increase in energy is accompanied by increase in the 
number of nodes in the eigenfunction, as illustrated in Fig. 3.8(b-d). In order to fully 
characterise them we need to assign two quantum numbers associated with them: 









(ne) 



220 4.74 728 9lB2 1226 1400'220 4.74 728 9B2 1236 14S0 

R (au) R (au) 


Figure 3.7: Probability density contours superimposed on the potential grid for eight 
resonance wave functions (a-h) at energies below Eth- The eigenenergies, En, are indicated 
in the box as E. 





E=0.9661 eV 


i E=11456 eV 


(c E=L25i6 eV 


(0) E=2.3650 eV 




Figure 3.8: ProbabiKty density contours superimposed on the potential grid for four res- 
onance wave functions (a-d) at energies above The eigenenergies, are indicated 
in the box by E. The number of nodes along the hyperspherical angle are: 6, 7, 8 and 18 
for a, b, c and d respectively. The closely lying resonant periodic orbits at 0.8861, 1.0649, 
1.2216 and 2.3470 eV «ire superimposed on the corresponding resonance eigenfunctions to 
illustrate the classical-quantal correspondence. 
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where ui corresponds to the hs^jerspherical radial (p) motion and 1/3 to the angular (^) 
motion. None of Fig. 3.8(a-d) corresponds to an excited state along p. Therefore we assign 
states (0,6), (0,7), (0,8), and (0,18) to the eigenfunctions plotted in Fig, 3.8(a)-3.8(d) in 
that order. It is worth pointing out that with increase in energy the hyperspherical modes 
are centered at successively larger p values and they cilso exhibit degeneracy. Understand- 
ably the lack of symmetry in the HeH^ system prevents us from assigning parity to these 
eigenstates. 

One thing that became clear from the examination of these eigenfunctions is that some 
of them show a tendency to decay into the (He-t-Hj ) channel, some into the (HeH'''+H) 
channel and some into neither of the two. Fig. 3.9 shows a function which reveals oscil- 
lations characteristic of clcissical trajectories for large R values and Fig. 3.10 reveals a 
hyperspherical mode which lives for a long time. 

The break-up of the hj^erspherical mode into motion that is dominantly in the R or 
r channel is illustrated in Fig. 3.11(a) for an energy slightly below and into motion 
along the dissociative channel as well in Fig. 3.11(b) for an energy slightly greater than 
the E^. 

Except for the two true boimd states shown in Fig. 3.7(a) and 3.7(b) all other eigen- 
states discussed in this chapter are of quasibound nature and they have, finite lifetimes. 
The Lifetime varies with the state. For example. Fig. 3.12 shows a spectral peak fitted to 
a Lorentzian (Eq. (3.9)) lineshape function for the resonance eigenstate with En = 0.9661 
eV which gives a value (Eq. (3.10)) of 215 fe. Other resonances at 0.9600, 1.2516, 1.3540 
and 2.5336 eV were found to have r„ = 140 , 445, 55 and 366 fs respectively. The one at 
iJ„=2.4442 eV (Fig. 3.p) has a particularly long lifetime of 1.91 ps. Unfortimately there 
are no other quantal estimates of lifetimes available for collinear (He,H2 ) collisions for us 
to compare our results with. 
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R (au) 

Figttre 3.9: Probability density contours superimposed on the potential grid for the reso- 
nance eigenfunction corresponding to the eigenenergy = 1.3540 eV showing the oscillation, 
characteristic of classical trajectories. 





72 



R (au) 


Figure 3.11: Probability density contours of resonance eigenfunctions superimposed on the 
potential grid at energies (a) slightly below and (b) slightly above Eigenenergies are 
indicated in the box by E. 
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3.3.2 Bound States and Resonances in Collinear HeHD"'' and 
HeDH+ 

The power spectrum, as obtained from time evolution of WP No. 2 (Table 3.1) for collinear 
HeHD'*' is shown in Fig; 3.13 along with the autocorrelation function |C'(f)| as an inset. 
The recurrences in |C(f)| persist until about 3-4 ps. This is reinforced by the probability 
density plots (not shown) on the potential grid during time evolution, carried out for a 
total of 5.29 ps. As mentioned earlier, the wave function is absorbed by the NIP (Eq. 
(3.8)) applied at the grid edges, bringing the norm < > down to < 10“^ at the end 

of the time evolution. The peaks in the spectrum correspond to quasibound states and 
hence to the resonances of the system. Although some of the peaks appear broad due to 
the overlap with the neighbouring peaks, many others have well resolved structure. T his 
is to be contrasted with the finding on the HeHj system, for which the the spectrum wais 
much more complicated. For HeHD"'" the resonances are less in number and are farther 
spaced. 

The transition state spectrum, obtained from time evolving WP No.2 (Table 3.1) along 
with the 1 0(^)1 as an inset for HeDH"'" is plotted in Fig. 3.14. The recurrences in 
|C(f)| persist till the end of time evolution (5.29 ps), by which time the norm of the wave 
function has come down to 0.027 . This implies that the lifetime of the HeDH"'" complex 
is more than that of HeHD"''and HeH^ . It is worth pointing out that the spectrum in 
Fig. 3.14 is more complicated than that in Fig. 3.13 and that the number of quasibound 
states is larger and more closely spaced for HeDH"'" than for HeHD"'" and HeH 2 "'' (cf. Fig. 
3.4). Interestingly, the spectra in Fig. 3.13 and Fig. 3.14 share some common peaks, 
particularly in the low energy region. This is due to the fact that those peaks correspond 
to the reagent diatomic vibrationed thresholds which are the same for both systems. 

The spectra obtained from the time evolution of WP No. 1 (Table 3.1) for HeHD "'"and 
HeDH"'" are shown in Fig. 3.15(a) and 3.15(b) respectively. Clearly, there are two bound 
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ECeV) 

Figme 3. 13: Transitioa sb^te spectnim in the energy range 0-3 eV for collinear HeHD''' 
:oinputed using tke Initial Gaussian wave packet INo. 2. Corresponding autocorrelation 
iinction is shown in thes inset. 
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Figure 3.14: Same as in Fig. 3.13 for coUinear HeDH^, computed using the initial Gaussian 
wave packet No. 2. 
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Figure 3.15; Transition state spectrum in the energy range -0.5 - 0.0 computed with the 
initial Gaussian wave packet No. 1 for HeHD'*' (panel a) and for HeDH'*' (panel b). The 
intensity 1(E) is plotted in arbitrary units. 
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states for both systems at E = -0.1507 and -0.0664 eV. The spectra differ only in terms 
of the intensity of the peaks. It is worth adding that our investigation on collinear HeHj 
had also revealed only two eigensates at energy below zero, but at -0.1281 and -0.0328 eV. 

Probability density contours of twelve resonance eigenfunctions of HeHD"'' 

are shown in Fig. 3.16(a-l). The eigenfunctions in Fig. 3.16(a) and 3.16(b) correspond to 
the two eigenstates having energy below zero. They correspond to the ground vibrational 
{v = 0) state having no nodes and the first excited vibrational (y = 1) state having one 
node of the complex. The wave function in Fig. 3.16(c) having two nodes corresponds to 
the V = 2 state of the complex of energy 0.2616 eV implying that there are no other bound 
states for collinear HeHD^. It corresponds to up = vp ( i/p corresponds to the number 
of nodes along the product channel (HeH+ -t- D) and up corresponds to the number of 
nodes along the reactant channel(He + HD''')) progression. The wave functions in Fig. 
3.16(d-f), correspond to up = up — Z progression. In terms of {up,up) assignment ( the 
spectroscopic u^ assigned to the HeHf resonances corresponds to up + up and ui is zereo 
for these resonances), these are (3,1), (5,2) and (6,3) respectively. This is in contrast to the 
findings of Skodje et al [100], who obtained the up= up — 1 and up = up — 2 progressions 
for the D-f-H 2 resonances. They did not find any up = up progression of the resonances. 
The two wave functions in Fig. 3.16(g) and 3.16(i) correspond to up = up — A progression 
and their assignments are (7,3) and (8,4) respectively. The eigenfunction in Fig. 3.16(h), 
in contrast, has a slightly distorted structure and can be assigned (7,4). The resonances 
in Fig. 3.16(j) and 3.16(k) are (12,6) and (14,7) respectively which belong to up = up/2 
progression, as is the case with the resonance in Fig. 3.16(i). Finally, the resonance in 
Fig. 3.16(1) corresponds to the (15,8) assignment and is more reactant channel specific. 
Therefore it becomes clear firom the above discussion that resonances in HeHD''' do not 
follow any unique progression. There are many other resonances for the system, but they 
could not be given any clear assignment owing to their complicated structure. We must 
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add here that all the resonances discussed above exhibit the characteristics of hyperspher- 
ical modes. 

Id. Fig. 3.17(a-l) the probability density contours of twelve resonance eigenfunctions of 
HeDH'^'are presented. Fig. 3.17(a) and 3.17(b) are the eigenfunctions of two bound states 
having the same energy as those of HenD"^. Fig. 3.17(a) represents the ground vibrational 
state (v = 0) and 3.17(b) the first excited vibrational state {v = 1) of the HeDH+complex. 
Interestingly the eigenfunctions for HeDH+ reveal a local mode character. These can 
be labelled in terms of quantum numbers representing He-D stretching iyne.—D) and H- 
D'^ stretching {vh-d*-) respectively. In terms of the above specification, Fig. 3.17(a) and 
3.17(b) correspond to (0,0) and (0,1) states respectively. Similarly, 3.17(c) and 3.17(d) 
represent the (0,2) and (0,4) states, in that order. With increasing energy the nature of 
the eigenfunctions becomes more complicated and any clear assignment becomes difficult. 
For example, the eigenfunctions in Fig. 3.17(e) and 3.17(f) exhibit an inverted ‘U’ shaped 
structure and the one in Fig. 3.17(g) shows a similar ‘U’ shaped structure and also the 
characteristics of a hyperspherical mode. The latter suggests ai/p= ur — S progression and 
can be labelled as (5,2) in terms of {i'r,i'p} as was done for HeHD"^. The eigenfunction 
in Fig. 3.17(h) also can be assigned (5,2). The one in Fig. 3.17(i) has a large probability 
density build up in the product channel (HeD"*" + H), and is difficult to assign. The wave 
functions in Fig. 3.17(j) and 3.17(k) have a matted structure. We have noticed a similar 
behavior in the case of Hj resonances [102] which are discussed in chapter 4. However, the 
nodal pattern along the hyperspherical radius (p) is nonuniform, particularly in the case 
of Fig. 3.17(k). If we consider only the nodes along the hyperspherical angle {^) these two 
states can be assigned (11,7) and (13,8) respectively, in terms of {vR,up) nomenclature. 
The eigenfunction in Fig. 3.17(1) seems to be quantized almost exclusively along the tj) 
direction and it contains a total of 25 nodes and can be assigned as (15,10). Many more 
resonances are found to have complicated eigenfunctions and are not shown here. Actually 









Figxire 3.17; Same as in Fig. 4, for HeDH+ resonances 
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HeHD- 


HeDH+ 


(eV) 

■Tn (fe) 

En (eV) 

Tn (fs) 

0.2616 

383 

0.9149 

510 

0.8074 

208 

1.0352 

645 

1.0142 

647 

1.5149 

136 

1.3443 

52 

2.0675 

480 

1.5086 

380 

2.3108 

454 

1.8679 

22 

2.6385 

378 

2.2230 

256 



2.4546 

514 



2.6599 

60 




Table 3.4: Eigenenergy (£?„) and line-width lifetime (r„) of some of the HeHD^ and 
HeDH"*" resonances. 


the large density of states makes it nearly impossible to identify all the resonances sys- 
tematically. As was found for collinear HeHj [268] this could be an indication of the 
underlying classical chaos. 

The lifetimes, of many of the resonances have been estimated from the width of the 
spectral peaks fitted to the Lorentzian (Eq. (3.9)) and they are listed in Table 3.4. Some 
of the peaks could not be fitted to a Lorentzian because either they have a 8- function 
shape or they are not true resonances (but a superposition). 

3.3.3 Vibrational State-Selected Reaction Probabilities 

The reaction probabilities over a wide range of energies extending up to the three body 
dissociation threshold (2.78 eV) for the reaction 

He + HD+(u) HeH+(i/) -k D (R2) 

for the initial vibrational states v = 0—3 of HD''' are plotted in Fig. 3.18(a-d). In all ceises, 
the dynamical threshold coincides with the energetic threshold. The reaction probability 
increases sharply with energy and then levels off; it rises again and levels off once again. 



Figure 3.18: State-selected reaction probabilities obtained by time-energj’’ mapping of the 
reactive flux of the wave packet for different initial vibrational states (t?) of HD for the 
reaction He + HD+(u) HeH^ + D. The threshold energies of various channels are 
marked in the lower abscissa of panel a to illustrate the occurrence of sharp variations in 
near each channel threshold. 
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This process is repeated until If{E) approaches a limiting value. This results in an 
ascending staircase like structure for the curve, particularly for u = 0 and 2 (Fig. 

3.18(a) and 3.18(c)). Fach step in the staircase occurs near the onset of a vibrational (t/) 
threshold of HeH"*" as can be seen from the threshold values marked on the lower abscissa 
of the Fig. 3.18(a). The steps in P^[E) in Fig. 3.18(a) clearly show that the prominent 
resonances for reaction (R2) are of threshold origin [277]. In addition to these threshold 
resonances, there are other resonances in the transition state spectrum (Fig. 3.13) and 
they can be correlated with the undulations of smaller amplitude in the plateau regions 
of the probability curve. The if (£) curve for v = 1 (Fig. 3.18(b)) gives no indication 
of resonances except one near the v' = 1 threshold. The magnitude of the probability for 
V = 1 remains nearly constant over a wide energy range. A similar behavior is seen in 
Fig. 3.18(d) for u = 3. It may be added that the if (^EO values are slightly lower than the 
if (£) indicating a slight vibrational inhibition. There is a marked drop in FfiE) when 
compaired to Ff{E) near the reaction threshold but at higher energies, if (£) is clearly 
larger than Pf{E). 

Dynamics of heavy-light-heavy reactions has received considerable attention over the 
years [278]. Oscillatory behavior of if (F) and the reactive scattering resonances for such 
reactions in collinear geometries have been well documented in the literature. An expla- 
nation for a step-wise behavior in cumulative reaction probability for the reaction in three 
dimensions in terms of quantized transition state energy levels has been proposed by Chat- 
field et al [279]. To the best of our knowledge, this is the first time that a clear cut staircase 
like structure for P^{E) is being reported and identified with threshold resonances. 

The reaction probability curves for the reaction 

He + DH+(i;) HeD+(t/) + H (R3) 

are reported in Fig. 3.19(a-d) for v = 0-3. The if(£) curve in each panel exhibits a dense 
oscillation, suggesting the existence of a large number of very closely spaced quasibound 






88 


states at the transition state. This is in agreement with the densely packed structure of 
the transition state spectrum reported in Fig. 3.14. Although the oscillations in P^{E) are 
sumlar for all the v states, their magnitude differs considerably. The probability values for 
V = 0 of DH"*" (Fig. 3.14(a)) are close to zero. Snapshots of the WP during time evolution 
reveal that much of the input translational energy gets converted into the vibrational 
energy of the inelastically scattered DH+. The reaction probabilities for u = 1 (Fig. 
3.19(b)) are about 10 times larger than those for the u = 0, If(£^ values are about 4 
times larger than JP^(jE/) and values, in turn, aure about a factor of 2 larger than 

FfiE). Therefore it is clear that there is an overall vibrational enhancement of reaction 
(R3) in contrast to the vibrational inhibition noted earlier for (R2). 

Perhaps the most striking feature of the Ff{E) values for reaction (R3) is the large 
number of oscillations, suggesting "chaotic" behavior in dynamics. While the term "chaos" 
seems to be well understood in classical mechanics, its meaning and implication in quan- 
tum mechanics continues to be the subject of debate. Sathyamurthy and coworkers have 
reported [263] on the chaos and fractals and their possible relation to quantal reactive 
scattering resonances for the collinear reaction (Rl) . More recently we have tried to char- 
acterize the quantal chaos in this system by analysing the transition state (TS) spectra 
and the results are documented in chapter 5. While many of the peaks in the TS spectrum 
(Fig. 3.13) for collinear HeHD"^ could be interpreted in terms of threshold resonances for 
the reaction (Rl) no such explanation is readily available for the TS spectrum of collinear 
HeDH"^. There seems to be a direct link between the oscillations in Ff{E) and the large 
number of oscillation in the TS spectrum in Fig. 3.14. 

3.4 The Negative Imaginary Potential Anomaly 

In section 2.3 we have discussed the applicability of NIPs when used in conjunction with 
the numerical time evolution schemes. In this section we compute the if (F?) values for 
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Figure 3.20: State-selected reaction probabilities for the reaction He -I- (o = 0) -)> HeH+ 
+ H, (a) computed by adding NIPs to the real potential (solid curve) and (b) computed 
by using a masking function (solid curve). The TIQM results for the same (R«f. [269]) are 
plotted as dashed curves. 
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the reaction (Rl) both in presence as well as in absence of NIPs. We have used the FFT 
method for the spatial evolution of the WP and the SO method for its time evolution. The 
SO method revealed an exponential damping of the WP in presence of the NIP in section 
2.3. However the Hamiltonian, becoimng complex in presence of a NIP, does not preserve 
the time reversal symmetry of the TDSE and possibly results in an anomaly. 

In Fig. 3.20(a) we show the Fq{E) values for the reaction (Rl) (solid curve) computed 
in presence of the NIP (eq. (3.8)) superimposed with the time-independent quantal values 
of the same reported by Sakimoto and Onda [269]. We have used the <Po{t) corresponding to 
the extended Rydberg potentials of Hj, computed by the FGH method [241], in the initial 
wave function (Eq. (2.124)). The same calculation wjis repeated several times by varying 
the height amd width of the NIP. The F^{E) values show the artificial oscillation at higher 
energies near the onset of collision induced dissociation. The P^{E) values computed by 
using the masking function (Eq. (3.13)) for damping do not reveal such oscillations and 
are presented in Fig. 3.20(b) (solid curve). The latter Pq{E) values are also in agreement 
with those obtained by Sakimoto and Onda (dashed curve). The i^(E) values obtained by 
Balakrishnan and Sathyamurthy [167] using a FFT-SIL propagator also do not reveal any 
oscillation and are in excellent agreement with the latter calculation. Balakrishnan and 
Sathyamurthy have used a damping function derived firom the Neuhauser and Baer’s lineeir 
absorbing potential*(Eq. (3.6)) and multiplied it with the time evolved WP at the end of 
each time step without adding it directly to the real potentials at the grid points. The 
damping function derived in this way is real and exponential and preserves the Hermitian 
property of the Hamiltonian. 

The implication of the above results is that the artefact in the numerical calculations 
result, particularly at higher energies, when NIPs are added to the real potentials at the 
grid points. The artefact disappears when 6in exq>onential damping function resulting from 
NIPs = -iVo) or any other real damping function (e.g. the one in Eq. (3.13)) 
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is used. This way of damping can be achieved more easily than the use of NIPs. The 
prescription outlined by Taylor and coworkers [236] on the use of the NIPs in conjunction 
with the CP scheme follows the latter idea amd is free from anomaly since the Hamiltonian 
never becomes complex. 

3.5 Summary and Conclusion 

We have computed the transition state spectrum of all three transients, HeHj, HeHD+ and 
HeDH"^ by the spectral quantization method and the reaction probabilities (IfiE)) for 
(HejHD"^) and (He,DH''') collisions through a time-energy mapping of the reactive flux of 
the WP. We have found that the spectra for HeHj and HeDH'*' are dominated by a large 
number of TS resonances compared to that for HeHD~, for which the threshold resonances 
seem to dominate. There are two boimd states found for all three transients. The analysis 
of resonance eigenfunctions from their nodal progression reveals the local mode behavior 
at lower energies and the hyperspherical mode behavior at higer energies. A comparison 
of the transition state spectrum of the three transients reveals that the density of state is 
Itirgest for HeDH"^. The recurrences in the |C'(t)| persists till the end of the time evolution 
for HeDH^ implying the existence of long lived resonances for this system when com- 
pared with the other two. The classical resonance periodic orbits are superimposed with 
the quantal resonances of HeH^ to establish the classical-quanta! correspondence. The 
energy resolved reaction probabilities for reaction (R2) reveal a characteristic staircase like 
structure which can be related to the threshold resonances, particularly for t; = 0 and 2. 
In contrast, the for reactions (Rl) (computed by Balakrishnan and Sathyamurthy 

[167]) and (R3) are highly oscillatory in keeping with their densely packed transition state 
spectra. This also suggests the possible chaotic dynamical behavior of (R3) for which the 
classical trajectory dynamics would be worth investigating. 



Chapter 4 


Transition state Resonances and 
Dynamics of Collinear (H“, H 2 ) 
Collisions 

4.1 Introduction 

Reactions involving hydrogen atom and hydrogen molecule and their ionic counterparts 
are of fundamental interest to the molecular dynamicists as they represent some of the 
simplest prototype exchange reactions. Two of them, viz., H + H 2 and H"'' + H 2 have been 
studied extensively both experimentally and theoretically [106, 279, 280]. However, much 
less is known about the reaction involving H": 

H- + H 2 ^ H 2 + H- (R4) 

The available experimental results on the system [281-283] provide an impetus to theo- 
retical studies, particularly in view of the competing electron detachment channel (with 
a threshold energy of ~ 1.45 eV) which becomes dominant over the charge transfer 
rearrangement channel for which the cross section becomes appreciable at ~ 0.45 eV. In 
contrast to the (H‘'‘,H 2 ) system which shows large vibrational inelasticity the (H~,H 2 ) colli- 
sions are only weakly inelastic [284]. Therefore, a comparative quantum mechanical study 
of the vibrational inelasticity in (H''’,H 2 ) and (H“,H 2 ) collisions have been imdertaken on 
the available ground state ab initio PESs for both the systems [285]. 
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EarHer calculation of the potential-energy surface (PES) [281] for the H 3 system showed 
a quahtative similarity with the neutral H 3 . Xhis was i^nforced by a more recent ab 
initio calculation of the ground state PES [286], which revealed a saddle point at = 

ri-H =1-997 a.u. with a barrier height of 0.4648 eV, relative to asymptotically separated 
reactants. The PES exhibits a van der waals ToiniTniTm in the collinear geometry which has 
been reported [286] to be capable of supporting four bound vibrational levels. Belyaev et 
al [287] computed the state-resolved reaction probabilities for the collinear H~ + H 2 (D 2 ) 
rearrangement reaction by S matrix Kohn variational method on a diatomics-in-molecules 
(DIM) PES and they did not reveal the existence of any resonance although a number of 
transition state resonances (TSRs) have been predicted for H 3 [288, 97, 98]. 

Therefore we have undertaken a study of the possible existence of TSRs in collinear 
(H“, H 2 ) collisions, by following the temporal evolution of wave packets chosen in the 
transition state region on the recently published ab initio PES [286]. We find that there 
exist a number of TSRs, at energies above the reaction threshold. While some of them 
could be related to the onset of higher vibrational charmels for the product H 2 , some others 
could not be. By computing the corresponding eigenfunctions we have obtained an insight 
into the nature of these resonances. 

We have also computed Ff{E) values for the collinear reaction (R4) to further check 
the existence of TSRs. The P^{B) values for the collinear reaction 

H + H 2 (v) H 2 (t/) + H (R5) 

are also computed on the Liu-Siegbahn-Truhlar-Horowitz (LSTH) PES [289] in order to 
compare the dynamical behaviors of the two system. 
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4.2 Computational Details 

4.2.1 Transition State Spectrum 

Since the detailed, methodology has been, described in previous two chapters we present 
only the essential details here. We have used mass-scaled Jacobi coordinates (12, to solve 
the TDSE on a 256 x 256 grid with the origin at (1.5, 0-5) a.u. in (JJ, r) with AJi = 0,02 
a.u. and Ar = 0.016 a.u. 

An even parity GWP (Eq. (3.2)) in H~-H and H-H stretching coordLuatess, - is plaaced 
initially in the interaction region of the PES and time evolved like for a half-collision. 
The spatial propagation is carried out by the FFT technique and the temporal evalintlon is 
carried out by the SO method. The time autocorrelation function is computed at eaclt ime 
step and half-Fourier transformed after filtering through the normalized IHa.iiniiigw"iQdow 
function to generate the transition state spectrum. The energy values conespc»ndLiiig to 
the maximum of the intense peaks in the spectrum are identified as the eigen-vaLues of 
the quasiboimd states that lead to resonances. The eigenvalues are accnratiel^f esstionated 
by fitting the spectral peaks to a line shape function assuming a single lioe fit_ TThe 
eigenfunction (^„) at any given eigenenergy (E„) is computed by Fourier tr'anasfcamaing 
the time evolved wave function at that particular energy. Lifetimes (t,^) of soame of the 
quasi-bound states corresponding to the resonances are calculated hy estimating tie full 
width at half TnaviTnnnri (FWHM) (r„) of the spectral peaks fitted to a loxentzian- fti_nc*ion 

(Eq. (3.9)). 

The reflection of the WP from the grid edges at longer times is prevented l>yaabsorbing 
it near the edges with the use of NIPs (Eq. (3.6)). The height and the width. oEthieMlP 
is appropriately chosen using the inequality condition (Eq. (3.7)) so as to ebri;ai3i99% 
absorption of the WP. 

We have carried out seven sets of WP calculations for the system by varying tbe location 
and the width of the GWP in the interaction region of the PES. The grid paramet ers for 
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those are listed in. Table 4.1. The first four sets of those calculations were performed 
on a finer grid with AJK = 0.02 a.u. and Ar = 0.016 a.u. and the remaining three 
were performed on a relatively coarse grid with AR = Ar = 0.05 a.u. We set the width 
parameter (5) of the GWP 0.25 a.u. except for the second set for which S was taken to be 
0.20 a.u. The upper limit of the potential Vn^ set on the grid is indicated in Table 4.1. 
We have time evolved the WP for a total of 32768 steps with At for each step indicated 
in Table 4.1 in order to bring the autocorrelation function down to zero. 

4.2.2 Computation of Reaction Probabilities 

The energy resolved reaction probabilities lf{E) are computed for the collinear reactions 
(R4) and (R5) on a 256x256 grid with the origin at (1.5, 0.5) a.u. in {R,r) and AR = 
Ar = 0.05 a.u. The Starck-Meyer (SM) PES [286] is used for (R4) whereas the LSTH 
PES [289] is used for (R5). The width parameter of the GWP is chosen (between 0.25 - 
0.15 a.u.) sufficiently narrow in coordinate space so that in momentum space it is broad 
enough to contain most of the energy components in it. This also exploits the power of 
the time-dependent wave packet approach in that a single WP propagation gives rise to 
reaction probabilities at various energies. The Morse oscillator wave function is used for 
(i>^{r) corresponding to the vibrational state of H 2 with vibrational quantum number v. 
The GWP is initially centered at iZg = 10.0 a.u. which is sufficiently far away from the 
interaction region so that initially the separation of variable condition (Eq. (2.124)) is met. 
The dividing line in the product channel is taken at rj = 8.5 a.u. The upper potential on 
the grid is fixed at 4 eV which allows At = 0.1616 fe for a stable propagation of WP by the 
SO scheme. The WP is time evolved for 662 fe which is found to be sufficient to remove 
essentially all the WP from the interaction region. The components of the WP reaching 
grid edges are absorbed by using NIPs (Eq. (3.6)). 
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WP No. 



< E >(eV) 

At 


1 

1.907 

2.247 

0.7142 

0.1077 

5.0 

2 

2.220 

2.526 

1.0909 

0.1347 

4.0 

3 

2.825 

3.209 

2.0017 

0.1077 

5.0 

4 

3.111 

3.461 

2.4910 

0.1077 

5.0 

5 

4.000 

4.000 

3.6482 

0.1077 

5.0 

6 

4.300 

4.300 

4.0204 

0.1347 

5.5 

7 

5.000 

5.000 

4.6093 

0.1077 

5.5 


Table 4.1. Parameters for different choices of initial wave packets used in time evolution 
in the spectral quantization technique 

4.3 Results and Discussion 

4.3.1 Transition State Resonances 

We illustrate the temporal evolution of the wave packet on the PES by plotting the prob- 
ability density contours, superimposed on the potential grid at various time intervals in 
Fig. 4.1. The initial location of the WP corresponding to WP no. 2 is shown in Fig. 4.1a. 
In 2.69 fe the WP spreads along the asymmetric stretch coordinate of H3 (Fig. 4.1(b)) 
and in 4.04 and 5.39 is it acquires a highly stretched configuration along the eisymmetric 
stretch coordinate as can be seen from Fig. 4.1(c) and 4.1(d) respectively. The WP starts 
developing structures in 8.02 fs (which in turn shows recurrences in the |C'(t)) plot) and 
shows progression along the symmetric stretch line (Fig. 4.1(e)). At 410 fe all of the WP 
is found to be removed from the transition state region and this time is sufficient to lead 
to a converged spectrum. However we have time evolved the WP further for considerably 
longer time (2.2 ps) to ensure a finer resolution in E. 

The transition state spectra along with the autocorrelation functions as insets, corre- 
sponding to WP No. 2 and 3 are shown in Fig. 4.2. It is clear from the insets of the 
figures that the |C(t)| at the higher energy decays at a much faster rate than that at a 
lower energy and the number of resonances axe far more in the higher energy region (Fig. 
4.2(b)). The amplitude of the oscillation in the autocorrelation function decays to zero 
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within 0.4 

ps in Fig. 4.2(a) and 0.2 ps in Fig. 4.2(b) which indicates the short-time nature of the 
(H JH2) collision dynamics. The intensity of the spectral peak at a particular energy is a 
strong reflection of the amplitude of that state in the initial basis set, which can be altered 
by varying the location and the width of the initial WP. Due to this reason some lower 
intensity peaks at one < > of the WP become the most intense peaks at some other 

<E>. 

The decay of the [(7(^)1 along with the transition state spectrmn in the form of a 
^og\Ql{E) plot as inset, resulting from the time evolution of WP No. 5 is presented in Fig. 
4.3. This WP is initially centered along the symmetric stretch line — tjj-h = 4-0 

a.u.) of the transient and has the average energy 3.6482 eV. The most striking featme 
of the spectrum in Fig. 4.3 is that it reveals the existence of a large number of resonances 
for the system in the higher energy (2 - 4.74 eV) range. 

A list of computed resonance energies is given in Table 4.2. We see that near each 
threshold of H2(u) there is one resonance [277] and besides that, there are other resonances 
too. These eigenvalues are quite different from those computed by Sadeghi and Skodje 
[97, 98] for collinear H3 which made us conclude that the resonances in collinear are 
significantly different from those of the neutral analog. 

Our conclusion is supported by the analysis of some of the eigenfunctions and their 
lifetimes. Eigenftmctions of four such resonances, shown in Fig. 4.4(a-d) in the form of 
probability density contours reveal a progression of nodes along hyperspherical radius (p) 
as well as hyperspherical angle {4>) with increase in energy. We have tried to classify 
these states in terms of vi (radial) and 1/3 (angular) quantum numbers as was done in the 
previous chapter. In addition we also label the states in terms of local mode quantum 
numbers (n, m±) by utilising its connection with (i/i, 1/3) [274, 290, 291]. 

The eigenfunction at 0.7767 eV (Fig. 4.4(a)) has essentially the local mode character. 
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Figure 4.1: Time evolutioa of the (initial) Gaussian wave packet (<E> — 1.0909 eV and 8 
= 0.2 a.u.) in the form of contours of l^p, superimposed on the potential grid at various 
times (0, 2.69, 4.04, 5.38, 8.08 fe and 2.2 ps) are shown in a, b, c, d, e and f respectively. 
The norm (< ^1’$' >) of the WP retained on the grid are: 1.0, 0.999, 0.995, 0.981, 0.942, 
and ~ 10“® in panels (a)-(f), in that order. 
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Figure 4.2: Transition state spectrum in the energy range 0-3 eV obtained from two 
different initial locations of the WP. The average energy of the WP is indicated inside 
each box. Intensity, 1(E) is plotted in arbitrary units. Insets show the decay of the 

corresporidiag autocorrelation functions. 
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This is a (0,2) state in terms of {u^, 1/3) specification and (1,1+) in terms of (n,m±) 
specification. The + sign indicates that this is an even parity state. The other states in 
Fig. 4.4(b-d) are dominantly hyperspherical modes. The state at 1.3796 eV (Fig. 4.4(b)) 
is a (3,6) state in terms of the former and a (3,6+) state in terms of the latter specification. 
Accordingly the 1.7160 eV state (Fig. 4.4(c)) is (4,8) (and (4,8+)) and 2.8664 eV state 
(Fig. 4.4(d)) is (7,14) (and (7,14+)) respectively. In terms of the {vi,up) specification, 
where u = 1/3/2, and p=parity, the states in Fig. 4.4(b),(c) and (d) are (3,3g), (4,4g) and 
respectively. The PES of collinear H3 is symmetric so that, we have parity choices. 
Since the states reported in this study have originated from even parity ’WPs, they all end 
up having even parity. However, the apparent asymmetry in the eigenfunctions in Fig. 4.4 
arises from the asymmetry of the {R, r) coordinates. 

The lifetimes of the states in Fig. 4.4(a-d) as estimated from the widths of the spectral 
peaks fitted to Lorentzians are 38, 30, 18 and 6 fe respectively, reiterating the short time 
nature of the dynamics. 

Probability density contours of twelve high energy resonance eigenfunctions of the 
collinear on the extended grid are shown in Fig. 4.5(a-j). It is obvious from the 
nature of these eigenfunctions that with increasing energy the corresponding quasiboimd 
states becomes excited along the asymmetric stretch coordinate. This is indicated by 
the disappearance of the nodal progression along the symmetric stretch (p) coordinate and 
increasing number of nodes along the asymmetric stretch (^) coordinate. The eigenfunction 
in Fig. 4.5(a) is a (8,6) state in terms of {ui,vz) assignment and the other eigenfunctions 
in Fig. 4.5(b) - 4.5(j) are (9,18), (7,18), (7,20), (3,22), (2,22), (3,24), (1,24), (1,26), and 
(0,28) respectively. Interestingly, the nodal progression of the resonance eigenfunctions 
resembles that of the H3 resonance eigenfunctions at higher energies. 
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V 

E„(eV) 

Et^(eV) 

E„(ir 3 )(eV)“ 

1 

0.7238 

0.7397 

0.7767 

0.782 

0.886 

2 

1.0409 

1.1212 

1.3796 

1.269 

1.300 

3 

1.4832 

1.5118 

1.7160 

1.728 

1.713 

4 

1.7771 

1.8384 

1.8987 

2.0591 

2.159 

2.108 

5 

2.1916 

2.4013 

2.5348 

2.5981 

2.563 

2.496 

6 

2.6255 

2.7094 

2.8664 

2.938 

2.868 

7 

3.0457 

3.0526 

3.1762 

3.2077 

3.2424 

3.283 

3.214 

8 

3.3533 

3.4048 

3.5509 

3.6417 

3.598 

3.532 

9 

3.7757 

3.9079 

3.880 

3.829 

10 

4.0172 

4.1542 

4.129 

4.085 

11 

4.2450 

4.3834 

4.341 

4.310 

12 

4.4145 

4.4654 

4.514 

4.492 

13 

4.7432 

4.643 

4.633 

Dissociation 


4.747 



Table 4 . 2 : Eigenvalues (£„) of the collinear H 3 resonances on the Starck-Meyer poten- 
tial-energy surface [286]. Ej,'* values are the energetic thresholds for the various H~ - 1 - 
H 2 (u) channels. EnCiTs) values are the collinear H 3 resonance energies, for comparison. 
“Reproduced from Ref. [98]. 
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4.3.2 Vibrational State-Selected Reaction Probabilities 

The P^{E) values for (R4) for v = 0 state of thus obtained are plotted in Fig. 4.6 
(as a solid line). Clearly, there is a large number of oscillations in Fq(E), far in excess of 
what one would have expected, based on the results of Belyaev et al [287] on a DIM PES 
and also based on the well known results [155, 292] for colUnear reaction (R5) in the same 
energy range. In order to ensure that we have not obtained any spurious results because of 
the choice of {R, r) grid or the computational scheme, we computed the lf(E) values for 
reaction (R5) under identical conditions on the LSTH PES [289] and found them (included 
as dashed lines in Fig. 4.6) to be in agreement with the values published in the literature 
[155, 292]. Thus it becomes clear that indeed there is a large number of reactive scattering 
resonances for collinear (H“, H 2 ) in contrast to a much smaller number for (H, H 2 ) despite 
comparable E{, values for the two reactions. We must add that for both reactions lf{E} 
becomes appreciable at ~ 0.45 eV, rises to near unity and eventually becomes nearly zero 
around 3.0 eV. 

The most prominent resonance observed for (R5) is at 0.882 eV, near the threshold 
of the u' = 1 channel of H 2 (0.782 eV). The first substantial dip in Po{R} occurs at a 
higher energy (1.1 eV) in case of (R5) and is less sharp. Our calculations [102] using the 
autocorrelation function approach had predicted the resonance at 1.12 eV. Preliminary 
investigations indicate that this is a Feshbach (compound state) resonance arising firom a 
quasibound state supported by the vibrationally adiabatic potential (VAP) for u = 2 in 
hyperspherical coordinates [293, 274, 255]. Fig. 4.6 also reveals that at higher energies, 
near each (vibrational) channel threshold, there is a resonance as anticipated by Friedman 
and Truhlar [277]. Some of the congestion near each such threshold resonance could be 
attributed to the Feshbach resonances arising from quasibound states supported by the 
higher VAPs. Still, there are some more oscillations in lf{E) remaining to be accounted 

for and their origin is presently imder investigation. 
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E(eV) 


Figure 4.6: Energy resolved reaction probability, P^, as a function of energy , E, for col^ear 
H- + H2 (v= 0) reaction on SM PES (soUd Une). The P«(E) values of collinear H + H2(v=0) 
reaction on LSTH PES are shown by a dashed line to highlight the dynamical differences 
between the two systems. The arrow along the abscissa indicates the threshold for the 

electron detachment channel. 
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The SM PES [286] has a shallow well of depth ~ 0.05 eV at {R,r) = (5.75,1.52) a.u.. 
But, this is unlikely to be the source of the large number of resonances found for (R.4). As 
mentioned earlier, there is an electron detachment channel opening up at ~ 1.45 eV and 
one can anticipate that some of the rearrangement channel resonances would be quenched 
at higher energies. Wie have tried to include the detachment channel by using a NIP along 
the deta chm ent seam in the PES. The magnitude of is highly sensitive to the choice 

of NIP and a more careful investigation is necessary. 

We must add here that we have checked for convergence in J^(E) by varying the height 
and the width of the NIP and the results have been further verified by obtaining them 
using the masking function (Eq. (3.13)) in place of the NIP. 

4.4 Discussion and Summary 

We have shown the existence of a large munber of resonances for collinear reaction (R4) 
on SM PES [286] in contrast to their absence on the DIM PES [287] and for reaction (R5) 
on the LSTH PES [289]. 

Qualitatively the ab initio PES for collinear Hg looks very similar to that of H3. 
Quantitatively, there are slight differences in the barrier height (^^'(.(Ha) = 0.418 eV on the 
DMBE PES [294] and 0.4249 eV on the LSTH PES [289]; = 0.465 eV [286] and 

in the barrier location (r^_;y = = 1.76 a.u. for H3 [294] and = ’"^-^=1-99 

a.u. for H3 [286]). The transition state for Hg being slightly more extended than for H3 
is understandable from the lower stability of the former, arising from the added electron 
in the antibonding orbital. This would also explain the tendency for electron detachment 
in Hg . The autocorrelation function for Hg decays faster than that of the neutral analog 
presumably due to the slightly increased barrier at an extended configuration for the 
former. As was pointed out earlier in the text all the resonances in collinear H3 could be 
identified as of (vibrational) threshold origin [277] but for Hg', there are also resonances of 
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different origin. In addition, it should be pointed out that all the resonances in H 3 axe of 
the asymmetric stretch type [97, 98] while no such generalization can be made in the case 
of H 3 . However, the asymmetric stretch type resonances dominate for the latter at higher 
energies. A. classical mechamcai study might lead to a clear understanding of the origin of 
some of the congestion in the curve as well as in the transition spectrum of the 

collinear (H'jHo) collisions. Such a study will be undertaken in future. 

The barrier for the exchange reaction is the lowest for the collinear approach in (H~,H 2 ) 
as well as (H,H 2 ). Therefore, the resonances reported herein could be considered repre- 
sentative of (H“,H 2 ) dynamics in three dimensions and the dynamical differences found 
here for the collinear geometry could be indicative of the differences in the dynamics in 
three dimensions for reaction (R4). Although for exchange becomes larger for lateral 
{C 2 v) approach, the threshold for electron detachment is the lowest in that geometry {Elf = 
1.2, 1.32 eV for the C 2 V and the collinear approach respectively). Therefore collisions along 
non-collinear geometries can be expected to show a larger tendency to electron detachment 
in (H",H 2 ). a three dimensional study of reaction (Rl) is presently in progress. 



Chapter 5 


Quantum Chaos in Collinear 
(He,H 2 ) Collisions 

5.1 Introduction 

Quantum chaotic behavior - the distinctive quantal feature of systems that show dynamical 
chaos in the classical limit — has been shown to be significant in a variety of systems 
of chemical interest [295, 296]. The usual methods of analysis are statistical, wherein 
the eigenvalue spectra of systems are examined for short- and long- range correlations, 
following similar analytical tools derived from random matrix theory (RMT) [297-303] . 

Satisfactory and reliable application of these techniques typically requires a large set 
of data (energy levels), which is often not accessible experimentally. Apart from numerous 
applications to model Hamiltonian systems with two or more degrees of freedom, such 
analysis has been applied to a number of systems of practical interest, including nuclear 
spectra, [297-299] the complicated electronic spectra of polyatomic molecules [304, 305] 
and vibrational spectra of rare gas clusters [306]. 

Irregular behavior in scattering systems is well known from both classical and quantum 
points of view [307]. However, the application of quantum chaotic analysis to realistic 
scattering systems, as for example chemical reactions, have been few, primarily since such 
analysis becomes possible only if the underlying PES can support a large enough number 
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of bound or quasibound states. The ion-molecule reaction Rl satisfies this requirement 
to a great extent: it is well known that there is a large number of narrow resonances 
[249, 251-253,255,256,259-262,167,266] corresponding to a large number of closely spaced 
quasibound states in its vibrationally adiabatic PES. In chapter 3, we have analysed the 
nature of some of these resonances by computing the transition state spectrum of the 
system (Rl) through time-dependent wave packet calculations although a complete and 
detailed analysis of all the resonances proves to be both laborious and complicated, owing 
to the large density of states. 

The present work re-examines the transition state spectrum of the above reaction using 
the methods of analysis of quantum chaos. 

Short-range correlations in quantum spectra can be studied through the nearest neigh- 
bor spacing distribution (NNSD). For systems which are integrable in the classical limit, 
the spacings have a Poisson distribution, while for systems which are classically completely 
chaotic, the spacings follow the Wigner distribution [297, 298]. Long-range correlations 
are examined through the spectral rigidity, Az{L) , or the two-point variance, S"(L), 
which again, are characteristically distinctive for regular and chaotic systems. 

The classical dynamics of the He,H 2 system is known to show chaotic behavior at 
sufficiently low energies [263]. A careful analysis [264] of the Poincare surface of section 
of the system reveals that chaotic regions of the phase space are surroimded by regular 
regions and vice versa. The quasibound spectrum originating from such a mixed phase 
space consists of energy levels from both regulcir and chaotic regions intertwined with 
the result that measures such as NNSD or A 3 (L) or TP'{L) reveal neither a completely 
regular nor a completely irregular behavior, but an intermediate between the two. Systems 
that are intermediate- neither completely chaotic nor completely regular- pose a problem 
inasmuch as interpolation between the two limits of behavior (of the NNSD or A 3 (L) or 
S^(L)) is necessary. Since many systems of practical interest fall in this category, [308] 
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this is a problem of considerable importance, and several suggestions on characterising 
their behavior have been made. These include a straightforward interpolation based on 
the relative fractions of regular and chaotic parts of the classical phase space, as well eis 
phenomenological forms [299, 308]. 

There is a further source of difficulty in these analyses. Although the HeHj scattering 
system has a fairly high density of resonances, the number is still far from ideal (or suf- 
ficient), and the sparse data set results in large fiuctuations in the NNSD, AsiL) and 
This problem can be overcome partially by examining other indicators of chaotic 
quantum behavior. The time domain characteristics of regular and irregular spectra are 
analysed by computing 

p{t) = \{m\m)? ( 5 . 1 ) 

which is the survival probability of a given initial wave function $(0). $(i) is defined as 

N 

. (5.2) 

i=l 

Following Pechukas, [309] Wilkie and Brumer [310] showed that there is a distinct 
difference in ((P(f))), the survival probability averaged over initial states and Hamilto- 
nian^ for regular and irregular spectra. Irregular spectra give rise to a ‘correlation hole’ in 
{{P{t))), which is a graphic evidence of spectral rigidity or level repulsion. This is absent 
in {{P{t))) for regular spectra. For systems that are intermediate in behavior Alhassid and 
Whelan [311] derived expressions for {{P{t))) . These methods have been seen to be of con- 
siderable utility in a quantitative analysis [312] of experiments on wave chaos in microwave 
cavities, where the number of states is insufficient for a clear diagnosis based on the NNSD 
or A 3 (T) alone. 

In the present work, we also explore the correspondence between classical trajectory 
dynamics and quantum dynamics of wave packets. Feit and Fleck [276] have studied the 
evolution of coherent Gaussian wave packets (GWPs) with a variety of mean positions 
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and momenta under the influence of Henon-Heiles potential and have shown the transition 
from regular to chaotic behavior with increasing energy (a similar kind of observation was 
made earUer by Hutchinson and Wyatt [313] ). They have investigated the phase space 
trajectories obtained from the expectation values of coordinates and momenta, behavior of 
the survival probability, nature of the power spectrum and the volume of the phase space 
covered by the wave packet during its evolution. All these diagnoses were found to lead to 
a reasonable interpretation of the system behavior. 

We perform s imi lar analysis for the reaction Rl. Three initial points in the phase 
space are chosen from the Poincare surface of section, [264] corresponding to the initial 
conditions of quasiperiodic, periodic and chaotic trajectories. The evolution of coherent 
GWPs corresponding to these initial conditions shows poor correspondence with the clas- 
sical dynamics, in contrast with the results of Feit and Fleck [276] for bound states. The 
difficulty in the present ca^e, we feel, is that this is a scattering situation with a relatively 
high density of states. The finite width of the wave packet extends over both chaotic and 
regular regions and thus gives rise to mixed or intermediate behavior. As may be expected, 
though, the correspondence improves significantly in the ^ ^ 0 limit. 

The spectral analysis of the transition state quasibound levels for He-Hj are described 
in detail in Section 5.2. The dynamical evolution of wave packets with different initial 
locations in the classical phase space is studied in Section 5.3.. The restilts are discussed 
and summarized in Section 5.4. 

5,2 Analysis of Eigenvalue Spectra 

The generation of the eigenvalue spectrum for collinear He,H^ system by time-dependent 
wave packet calculations hsis been described in detail in chapter 3. A typical spectrum, 
obtained from the time evolution of an initial Gaussian wave packet with an average 
energy, {E) = 0.9924 eV, placed in the interaction region of the MTJS PES [254, 258] is 
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shown in Fig. 3.4, The spectrum contains approximately 170 quasibound states up to 
the limit of the three body HeH 2 dissociation threshold (2.78 eV). The complex nature of 
the spectrum is due to the mixing of diflFerent frequency components of a large number of 
classical orbits of the system, which is an indication of the chaotic nature of the underl 3 dng 
dynamics. 

Prior to statistical analysis it is essential to “unfold” the spectrum, i.e. make the 
average density of energy levels uniform over the entire energy range. A new set of energy 

levels 1 ^ 1 , J52, j is calculated from the old set ,-En} by dividing all the 

level spacings by the local average spacing [314, 315] 

EiJri Ei + {2k 4- 1) — — ^ = — 0)\j2 = i + k,{i + k < n) (5.3) 

where k is the number of consecutive spacings used to calculate the local average spacing. 

The short-range correlation between the levels is estimated from the NNSD P{s)^ where 
s is the spacing between two consecutive levels normalised with respect to average level 
spacing which is, after unfolding, 1. That is, P{s)ds = probabiUty of finding two consecutive 
energy levels at a distance s apart from each other, in which the lower one is at x and the 
upper one is in the interval x + s<E<x + s + ds with the condition 

Jp{s)ds =1 = JsP{s)ds. (5.4) 

For a system known to be classically chaotic, the spacing distribution shows ‘level repul- 
sion’ with P(s) 0 as s — > 0 whereas, for a classically integrable system there is ‘level 
clustering’, P{s) 1 as s — »• 0. For the former, the distribution curve, in the semiclassical 
limit, is given by the Wigner surmise [298, 316]: 

P(s) = yexp (5-^) 

while for the latter, the Poisson distribution [298, 317] 

P{s) = exp{-s) 


(5.6) 
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is applicable. The spectral rigidity, AsiL) , and the two - point variance, E^(Z) , 
calculated from the unfolded spectrum, are used to diagnose long-range correlations [298] 
between the eigenstates in an interval of energy, L. A^^L) actually measures the least- 
square deviation of the staircase function the number of energy levels below energy E 

from the best straight line at any interval. Following Dyson and Mehta [318], Az{2L) can 
be written as 


As (2i:) = 



Min{A,B) [n{E) 


-AE-B 



(5.7) 


where AE^-B represents the best straight line fit of the staircase in the interval [x-L, x-\-L]. 
In terms of unfolded eigenvalues { Ei, Ez, Ez , ..., J5„} lying in the interval [-L, L] A3(2X)is 
given by [319] 


Az{2L) = 






+ 


2L 


'y ^ — 2i "t" l)jEi 

^ 1=1 y 


( 6 . 8 ) 

where n is the actual number of levels within each ensemble of the unfolded spectrum. It 
is well established [296, 298, 320] that the Gaussian orthogonal ensemble (GOE) limit of 
Hamiltonian matrices provides a suitable standard to measure chaos. 

The chaotic limi t of A,.z{L) is identical to the GOE value of RMT: 


{L) = \logL - 0.00695 


(5.9) 


On the other hand, A 3 (L) for an integrable system follows the Poisson value 

(L) = ^ (5.10) 

We have also computed E^{L] [298, 299, 301, 321] of the number, n, of levels in the 
interval x to x + L : 

S'‘{L) = ([n(i,i)- < n(x,L) >]“) 


(5.11) 
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E“(i) measm-es the two-pomt fluctuation iu the spectrum and it gires the same loug- 
range correlation as A^iL) . The GOE value of ^^(X) is given by 

^Gosi^) = 2^2 (X) + — [5i(7rX)]^ - i5'i(7rX) (5-12j 

TT TT \ > 

with 


^ 2 (-^) = ^ [logi2TrL) + 7 + 1 - cos(27rX) - Ci{2irL)] + X 


^S'i(27rX)j 


(5.13) 


Quantal fluctuation patterns in the time domain are analysed by computing the survival 
probability using Eq. 5.1. In case of regular spectra the survival probability, averaged over 
initial states and Hamiltonian, ((i^t))) , starts from an initial value 1.0 and approaches the 
asymptotic value monotonically. In case of an irregular spectrum, however, {{P{t))) falls 
below the asymptotic value at short times (t = 2Trh/ < AE >) and then approaches 
the asymptotic limit at a later time, thus giving rise to a ‘correlation hole’. To compute 
((P(i))), the averaging over the initial states is done through [322] 


(W) = 


^■ + 2 


1 + 


37V 


E 


COS 




27r(X’n — Em) 


t 


< AE> 


(5.14) 


and the averaging over the Hamiltonian is done by averaging over n levels which are 
divided into several segments, each containing N levels. 

The NNSD, P{s), for four different transition state spectra are shown in Fig. 5.1(a- 
d). Each spectrum originates from different parts of the interaction region of the MTJS 
PES, and therefore differs in detail, in terms of some of the eigenenergies and intensities 
of the common peaks. In all cases, the distribution is neither Poisson nor Wigner, but is 
intermediate in behavior. 

In Fig. 5.2 (a-d) the spectral rigidity, A 3 (X) , is plotted as a function of X. Again, the 
behavior is intermediate between the GOE and Poisson limits, although the fluctuations 
are large, resulting from the relatively small number of energy levels in the system. The 
two-point variance, (which is not shown here), varies in a similar manner. 
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Figure 5.1; The nearest neighbor spacing distribution, P{s), as a function of the spacing 
between energy levels, s, showing short - range correlation, (a), (b), (c) and (d) are for 
170, 163, 165 and 131 quasibound states obtained from four initial Gaussian wave packets 
with (E)= 0.9924, 1.2159, 1.7694 and 2.2816 eV respectively. Poisson and Wigner indicate 
the distributions in the regular and the irregular limits respectively. 







117 



Figure 5.2: The spectral rigidity, A 3 (L) , as a function of the segment l^gth, !,, showing 
the long - range correlation, (a), (b), (c) and (d) corre^ond to the WPs referred to m 
Fig. 5.1. Poisson and GOE indicate the regular and the irregular limit respectively. 
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The behavior of the average survival probability for the four sets of energy levels is more 
revealing, and is shown in Fig. 5.3(a-d). In our calculation we use N — 10, which gives an 
asymptotic value of ((i^t))) = 0.25. As can be seen from Fig. 5.3, at short times ((P(jt))) falls 
below the asymptotic value and then only it approaches the asymptotic limit. Again, the 
somewhat large fluctuations of at long times around the asymptotic value is due 

to the inadequate number of energy levels in the system. However, the signature of level 
repulsion is distinctly present in all the four sets of resonances. Due to the mixed (partly 
integrable and partly non-integrable) nature of the phase space the ((F(t))) curves show a 
behavior that is intermediate between the regular and chaotic limi ts, (For illustration, the 
behavior of {{Pit))) in the regular limit, estimated from a set of 8000 random numbers, is 
included in Fig. 5.3(a), as shown in dots.) 

Since ail the measures discussed above indicate that the system has mixed dynamics, 
we can estimate the fraction of the regular and the chaotic regions in the phase space by 
fitting the experimental {i. e. derived from the present data) ({P(t))) curve to that predicted 
from an RMT analysis. For mixed systems, {{P{t))) can be written as [311] 

mm = {l + I [«(t) - 6,,(r)l} (6.15) 

where r = t/27r < p > {< p > is the average level density) and ^ = 1, 2 and 4 for Gaussian 
orthogonal, unitary and symplectic ensembles respectively. ?) 2 / 3 (t) is the two-level form 
factor [298] obtained by Fourier transform of the two-point cluster function, ^ 2 ^( 0 ;). The 
asterisk signifies the convolution 

An *fit) = Jdt' An {t')f {t - 1') (5.16) 


and 


A;, (t) = jr> 


sm{i:Nt) 

irt 


(5.17) 


The convolution with Ajv results from the use of finite subspace, N. If a fraction, j3. 



Figure 5.3: Variation of survival probability, ((P(t))) with time t (a), (b), (c) and (d) 
have the same meaning as in Fig. 5.1. {{Fit))) for a regular spectrum (generated 
8000 uniformly distributed random numbers) is shown by a dotted Ime 1 °- Fig. 5.3(a). The 
Rl^IT ((P(*))) values for intermediate are included in each panel by a dashed hne. 
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of the levels obey GOE statistics and 1-/3 obey Poisson statistics, bos{t) {0 < /3 < 1) is 
given by 

heit) = 0b2{t/0) + [(1 - /3)/N] 6{t) (5.18) 

where b 2 indicates the exact GOE value which is given from the two-level cluster function 
[297, 298] as 

b 2 (t) = l-2|i| -h |t|log(l-f 2|tl) 1*1 < 1, 

= -1 + |i| log((2li| + l)/(2|t! - 1)) \t\ > 1 (5.19) 

Eqs. 5.17, 5.18 and 5.19, when substituted in Eq. 5.15 give the desired expression for 
((P(t))) for the mixed case, which is superimposed in Figs. 5.3(a-d), as a dashed line. The 
agreement is reasonable for a value of /3 « 0.9, indicating that the effective chaotic region 
is large at these energies. We return to this point in Section 5.4. 


5.3 Wave Packet Dynamics 


The dynamical behavior of a wave packet in coordinate space is studied by solving the 
TDSE in mass-scaled Jacobi coordinates {R, r). The H am iltonian for collinear He, Ho is 
given by Eq. (3.1) The initial wave function is taken as a coherent state GWP, 


$(jR, r,t — 0) = Nexp 


{R-RaY iPiR 

itfl h 


exp 


(r-ro)^ iP°r 
+ 


(5.20) 


2ct^ h 

where N is the normalization constant and cr is the width parameter of the wave packet. 
The initial conditions are specified in terms of the average quantities 


<R>=Ilo; <Pr >= Pr 

<r>=ro; <Pr>=P° (5-21) 

As has been described in chapter 3, the spatial propagation of the wave packet is carried 
out by the FFT method [66] and the temporal evolution is carried out by the SO 
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Wave packet 


Initial Conditions 


(i!) 

■ (eV) 

Behaviour 

(eV) 

< R> 

< r > < Pr > 

< Pr > 

(a) 

3.5780 

1.5623 -0.2191 

0.7902 

0.2148 

Regular 

(b) 

4.3360 

1.1865 0.0 

0.0 

0.7091 

Regular 

(c) 

4.2170 

1.5623 -0.0626 

2.0106 

1.2263 

Chaotic 


Table 5.1. Parameters for different choices of initial wave functions used in wave packet 
dynamics. 


method [64]. The parameters chosen for the initial wave function are listed in Table 5.1. 
All calculations were done on a 256 x 256 spatial grid with the origin at (2.2, 0.4) a.u. in 
(jR, r) space. The grid spacings are 0.05 a.u. in both R and r. The upper limit of the 
potential on the gnd is set at 3 eV which allowed a time step of 0.19 fe in the propagation 
scheme. The width parameter of the GWP is taken as 0.25 a.u.. The survival probability 
(Eq. (5.1)) is computed at each time step and the power spectrum is generated through 
Eq. (2.116). The volume occupied by the wave packet in phase space during its evolution 
is computed from the uncertainty product as [276] 

V(t) = |< {Pr- < Pr >f >< {Pr- < Pr >f >< {R- < R >Y >< {r- < r >Y >}5.22) 

The time evolution of three wave packets in (iZ,r) space is shown in Figs. 5.4(a-c). 
Initial conditions are chosen from the Poincare surface of section [264] of the system Rl. 
Although the behavior of the classical trajectories originating from the initial conditions 
corresponding to WP (a) and (b) has been found to be regular, the wave packets, because 
of their non-zero width, extend to nearby nonintegrable regions of the phase space and 
spread. The spreading of a coherent GWP is dependent on the extent of anharmonicity of 
the region of the PES it is located initially. The energy of the wave packet in Fig. 5.4(c) 
is large compared to that in Fig. 5.4(a) and 5.4(b) and since it is located initially at a 
position corresponding to a chaotic trajectory, there is a large degree of spreading. 











Figure 5.4: Probability density contours (1$P) of the initial coherent Qaussian wave packet 
superimposed on the potential energy surface showing spreading oi the wave pacicec wicn 
time, (a), (b) and (c) are for the different. initial conditions given in IaD.ie 5... 
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The volume V{t) of the phase space occupied (calculated through Eq. (5.22)) by the wave 
packets during their time evolution is shown in Figs. 5.5(a-c). A large increase in V(t) 
for the wave packet (c) compared to (a) and (b) is in agreement with its high degree of 
spreading in phase space. 

The shape of P{t) for the three wave packets is shown in Figs. 5.6(a-c) respectively. An 
exponential decay at very short times is observed due to the dephasing of the initial wave 
packet. This phenomenon is well known in quantum mechanics and the dephasing time is 
given by where F is the energy spread of the wave packet. Classically this is the 

time needed for the wave packet to move a distance equal to its width (t^ = <t/u, where v is 
the wave packet velocity). For time t > Td, P(t) exhibits a large number of recurrences of 
different amplitudes which is composed of a large number of independent frequencies, and 
hence accounts for the spreading of wave packets. The wave packets (a) and (b) exhibit a 
moderate correlation till the end of the time evolution and hence reveal intermediate kind 
of behavior with a high degree of reguleirity. The wave packet (c) on the other hand does 
not show any correlation which is an indication of the chaotic behavior. 

We have examined the power spectrum and found that the wave packets (a) and (b) 
show a moderate degree of correlation between the successive peaks and there is a high 
degree of regularity hidden inside the overall irregularity. On the other hand the spectrum 
of (c) is purely uncorrelated and shows the contribution of a large number of states with 
a high degree of anharmonicity. 

We have also studied the behavior of the wave packet (b) in the classical limit ( 0 ). 
Classically the initial conditions of wave packet (b) correspond to a periodic orbit. This 
wave packet does not show any spreading as h ^ 0. The behavior of P{t) also shows a 
high degree of correlation and the corresponding power spectrum is regular with constant 
spacings between successive peeiks. 
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5.4 Discussion and Summary 

The signatures of quantum chaos in the scattering dynamics of collinear He,Hf can be 
seen when the transition state spectra are analysed using commonly employed statistical 
measures such as the nearest neighbor spacing distribution (Fig. 5.1), the spectral rigidity 
(Fig. 5.2), and the averaged survival probability (Fig. 5.3). Due to the rather small 
number of states present in the system, there are large statistical fluctuations, but in all 
cases it is clear that the system shows level repulsion as well as evidence of being of mixed 
typ®: neither regular nor completely irregular. The extent of chaos in the system 

is estimated to be quite large. By fltting ((P(t))) to the RMT, we And that reasonable 
agreement obtains for /5 = 0.9, corresponding to 90% of the states being irregular. This 
result is not surprising since, as we have discussed earlier, each spectrum originates from 
a different part of the interaction region of the PES, and the levels are contained within a 
small section of the phase space. To have an idea of the overall irregularity of the phase 
space, we did a sepairate analysis with the whole spectrum obtained by combining the four 
spectra: {{Pit))) for this case is shown in Fig. 5.7(a), and is best fit by « 0.6, suggesting 
that approximately 60% of the quantal resonances in R1 are irregular in nature. 

In an earlier study, [263, 264] Sathysimurthy and coworkers have investigated the nature 
and extent of classical chaos in Rl. The fractional number of trajectories exhibiting chaotic 
behavior decreased with increase in energy ais shown in Fig. 10 of Ref. [263]. It is not 
possible to draw a one-to-one correspondence between classical chaos and its signature 
in the quantum dynamics as the former pertains to a particular energy while the latter 
involves a distribution of energies. 

Complementary results were given in Section 5.3, where we studied the correspondence 
between the classical trajectory dynamics and quantum wave packet dynamics. It has been 
shown by Bixon and Jortner [323] that such a correspondence is strictly valid only when a 
coherent GWP evolves under the action of a harmonic potential. In such situations motion 
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Figure 5.7: Variation of survival probability, {{P{t))) with time t (a) ioi the overall spec- 
trum of the phase space, constructed by combining individual spectra, and (b) for the 
spectrum obtained from the WP (b) illustrated in section 5.3. The KMT ((jR(t))) values for 
13 = 0.6 (for (a)) and for /3 = 0.7 (for (b)) are superimposed and indicated by the dashed 
lines. 
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of the wave packet is described by the first moments < g > and < p > and the center of 
the wave packet moves along the classical trajectory. However, appearance of terms first 
order in second moments induces spreading of the wave packet, and due to this the center 
of the wave packet moves away from the classical trajectory and the correspondence breaks 
down. 

The total voliime, V{t ) , of the phase space the wave packet occupies during its evolution 
is a quantitative measure of the extent of its spreading. It can be shown that the entropy, 
5, derived [324] from the partition function varies linearly with nfcglnVl Therefore, an 
increase in V gives rise to an increase in S, Our results in Fig. 5.5 show that the entropy 
associated with the evolution of the wave packet (c) is larger than that for wave packets 
(b) and (c). Hence, the behavior of the wave packet (c) can be reasonably classified as 
irregular and that of wave packets (a) and (b) as relatively regular. 

The change in the initial state population with time which is termed as the survival 
probability has been extensively used to check correlation. The semiclassical limit of P{t), 
as discussed in Section 5.2, is the ensemble averaged value. However, P(t) computed firom 
the wave packet itself contains all the recurrences and shows the following features: 

(i) a wave packet initiated at a fixed point corresponding to a stable periodic orbit shows 
a periodic variation of P{t) with constant amplitude of oscillation which denotes a high 
degree of correlation; 

(ii) for a wave packet initiated at a point corresponding to a quasiperiodic trajectory, P{t) 
falls sharply below 1.0 and then oscillates with reduced amplitude, signifying moderate 
correlation; 

(iii) for a wave packet initiated at a point corresponding to a chaotic trajectory, P(t) goes 
to zero sharply as the center of the wave packet moves away due to spreading and shows 
no correlation. 

The behavior of P{t) in Figs. 5.6(a) and 5.6(b) conforms to (u) above and thus the 
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corresponding motion can be classified as quasiperiodic, while the behavior in Fig. 5.6(c) 
conforms to the point (iii) above and hence the corresponding motion is chaotic. The 
behavior of the power spectrum is complicated in all the cases. However, we find that the 
power spectrum of wave packets (a) and (b) shows a high degree of regularity and that of 
(c) a large irregularity arising from the mixing of the different frequency components. 

RMT analysis of the power spectra obtained from (a), (b) and (c) does not reveal any 
distinct difference in the behavior of these spectra in the statistical limit. For example, 
the ensemble averaged survival probability, ((P(t))) , computed by RMT gives (3 = 0.7 for 
all of them. ((P(i))) generated from the power spectrum of WP (a) (solid line) along 
with ((P(i))) computed by RMT (dashed line) is shown in Fig. 5.7(b). Since the three 
wave packets differ marginally in terms of their location on the PES and all of them have 
the same width parameter, the initial energy components present in their basis set are 
essentially the same. The only difference exists in the weights of the energy components. 
This difference is reflected in the intensity of the peaks in the spectra. Since the statistical 
analyses do not depend upon the intensities but on the location of the resonances, all three 
spectra show similar statistical behavior. 

Although the quantinn wave packet dynamics reproduces the general classical behavior 
successfully, quantitative correspondence is lacking. The description can nevar be Complete - 
as- the uncertainty principle is at the heart of quantum dynamics and quantum interference 
effects have no classical analogue. Investigation of the behavior of the wave packet (b) in 
the classical liroit ( /i -> 0 ), shows that spreading is reduced in this limit and ‘scars’ in 
the wave packet coalesce around a point [325]. The power spectrum of P{t) shows peaks 
at constant spacings and very little mixing of states with different frequencies, verifying 
that the quantal-classical correspondence becomes pronounced in this limit, as expected. 

In summary, we would like to point out that quantal calculations [269, 167, 101] have 
revealed a large number of resonances for R1 and classical trajectory calculations have 
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shown a number of trajectories to be irregular [263, 264]. In the present study we have 
shown, using various measures as NNSD, A 3 (L)and J:\L) that the quantum dynamics in 
the system can. also be characterized as partly regular and partly irregular, thus establishing 
a direct hnk between classical chaos and quantal resonances in t his important, prototypical 
reactive scattering problem. 



Chapter 6 


Transition State Resonances in Three 
Dimensional (He,H^) Collisions 

6.1 Introduction 

A number of studies [326] have shown that reactive scattering resonances are important 
in three dimensional collisions. The peaks in the photodetachment spectra obtained by 
Neumark and coworkers [44-46,109] for a number of heavy-light-heavy (H-L-H) systems 
have been reproduced by the theoretical calculations of Schatz and coworkers [111,112,114- 
116] . They have been related to the topological features of the corresponding PESs and 
to the oscillations in the Ff{E) curves for the H -f LH collisions. Kuppermann eind 
coworkers [327] have computed the resonance eigenenergies and lifetimes for the 3D H 
+ H 2 (J = 0, 1) collisions on the PK2 as well as LSTH PESs. Based on the results of 
their computation they have anticipated that the 3D resonance energy can be obtained 
by adding the separable bend energy to the corresponding collinear resonance energy. 
Recently, Skodje et al [99] have computed the resonance eigenenergies and lifetimes for the 
3D H + H 2 ( J = 0) collisions on the DMBE PES by extending the SQM in three dimensions 
and found them to be in good agreement with the results of Kuppe rman n and coworkers 
[327]. They have reported 33 resonances in 3D [99], in contrast to 13 for the collinear 
geometry [97, 98] of the system. They have assigned quantum numbers to the 
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eigeafunctions. In addition they have also mentioned that many more resonances exist for 
the system in 3D, but they could not be assigned clearly because of their complicated nodal 
pattern and anticipated that the inclusion of the coupling with the higher excited state 
in the calculation is necessary at higher energies. It also became clear that the colUnear 
resonances exhibit progression along the bending coordinate in three dimensions, and the 
number is far in excess of that for the collinear geometry for the simplest prototypical 
exchange reaction. 

A large number of narrow resonances have been shown to exist in collinear (He,Ht) 
collisions (chapter 3) and quantum chaos has been inferred (chapter 5). The existence 
of resonances in 3D He+Hj‘(J = 0) collisions has been shown by Kress et al [256] and 
Lepetit and Launay [259] on the MTJS PES. They reported resonances having lifetimes 
up to ~ 3 ps. Results of Zhang et al [260] and Balcikrishnan and Sathyamurthy [265] for 
3D J = 0 collisions also suggested the existence of numerous resonances on the same PES. 
Recently, Mandelshtam et al [266] have correlated the quantum mechanical resonances to 
the classical chaotic scattering through an analysis of the periodic orbits. We have tried 
to investigate the resonances and their cheiracteristics in 3D He+Hf (J = 0) collisions by 
extending the SQM in 3D. While the methodology is described in section 6.2, preliminary 
results are presented in section 6.3 and a succinct summary of our findings follow in section 
6.4. 


6.2 Methodology 


The nuclear Hamiltonian for 3D A+BC ( J = 0) collisions in mass scaled Jacobi coordinates 


is given by [225]: 

H 


~2 ;2 


2/i 
2/i 


2/i 

BE? ^ 


1 d 
2Jsin7^ 



+ V(i!,r, 7 ) 


( 6 . 1 ) 
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where 7 is the angle between R and £7 = 0 corresponding to the approach of A to 
the B end of the BC molecule, /x is the three body reduced mass and I = p.{R? + is 
the moment of inertia and j is the angular momentum operator. The associated voliune 
element is dr = sin'^dRdrd'y. 

The initial wave function is taken as an even parity GWP, r, 7 , t = 0), in terms of 
the bond coordinates of the triatomic complex ABC: 


$(i?. r, 7 , f = 0) = Nexp 


{^AB - r%y (rBc - ' 




where N" is the normalization constant. The AB and BC intemuclear distances and 
rgc Stud the bending angle &abc are expressed in terms of (R, r, 7 ). < 745 , crgQ and xabc are 
the width parameter of the GWP along the respective coordinates. 

The TDSE is solved in (i?, r, 7 ) coordinates. The action of the radial and angular 
kinetic energy operators on the WP is accomplished through the FFT method [ 66 ] the DVR 
method [68.228-230] respectively. The sampling points along the 7 coordinate are chosen 
corresponding to the Gauss-Legendre quadrature points [211]. The grid representation of 
the initial WP at the point (E/,r^, 7 n) is given by: 


= ^(Rl, rm, 7n, t = 0) = Jn, t = 0) 


(6.3) 


where w„ are the Gauss-Legendre weights at 7 n- The temporal evolution of the WP is 
carried out by the SO method as follows [328]: 

+ 0[(A4)=) (6,4) 


In numerical calculations this is done as follows : 
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(i) a 2D forward FFT takes the WP from the {R,r) space to the corresponding 

(is, fe.) space where it is multiplied by the exponential containing the diagonal of 

the kinetic energy operator and then brought back to the (R, r) space through an inverse 
2D FFT; 

(ii) the results from the above step are then transformed to the an gular momentum (j) 
space (FBR) through the DVR-FBR transformation matrix 3^ where the angular kinetic 
energy operator is diagonal, multiplied by its diagonal values contained in the exponent 
and brought back to the discrete angle space (DVR) through the inverse transformation 
matrix 

(iii) since the potential energy operator is diagonal in the discrete space, the result from 
step (ii) is simply multiplied by the values of the exponential containing the potential 
V{Ri,rm-i'Yn) at the discrete grid points (1, m,n); 

(iv) steps (ii) and (i) are repeated on the result. 

In practical implementation the time axis is sliced into N steps of length At each as 
mentioned in Sec. 2.4.1, and the entire procedure (i) - (iii) is repeated for each At. The 
DVR-FBR transformation matrix in the Gauss-Legendre quadrature DVR is defined by: 

T;n = v^fli;(cos7„) (6.5) 

where (l>j{cos'y) are the orthonormal FBR basis functions (normalised Legendre polynomi- 
als) which are written in terms of the usual unnormalized Legendre polynomials, P;(cos7): 

(j)j{cos 7 ) = yS^Pj(cos 7 ) (6-6) 

The maximum number of angular basis functions used to construct Ty„ is equal to (iST.^— 1), 

where is the number of grid points along 7. 

As the WP different time steps is computed, the autocrrelation function is 

also calculated (Eq. (2.115)) and finally Fourier transformed to the energy domain (Eq. 
(2.116)) to generate the power spectrum. The peaks in the power spectrum are analysed 
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in the same way as described m chapter 3, to obtain the eigenvalues and the lifetimes of 
the resonances. The corresponding eigenfunctions are then computed through Eq. (2.122) 
and analysed to characterize the resonances in terms of {uuU 2 ,i'z) quantum numbers. 

In order to compute the He + Hj(J = 0) resonances on the MTJS PES we constructed 
a (128 X 128 x 16) grid in (i2,r,7) space with the origin at (2.2 a.u., 0.4 a.u., 0.145 rad) 
and AR = Ar — 0.07 a.u., A7 ~ 0.19 rad. The sampling points along 7 correspond 
to the 16 point Gauss-Legendre quadrature. The initial GWP is centered at = 

3.20 a.u., T^ff+ — 5.3 a.u. , ~ 2-99 rad) and the width parameters are taken as 

<^HeH ~ — 0.25 a.u. and = 0.15 rad. The angulair width is chosen such that 

appreciable bend excitation can occur in the complex HeHt. We set At = 0.1616 fe and 
the WP is time evolved for a total of 5.29 ps which gives a hi gh enough resolution (=7.8 x 
10'"'^) in the transition state spectrum. The components of the WP reaching grid edges at 
a later time are damped out by using the masking function (Eq. (3.13)) which is activated 
at Rj = 9.69 a.u. and r/ =7.89 a.u. and covers 16 % of the grid points in (i?,r). 

6.3 Results and Discussion 

The snapshots of the time evolved WP corresponding to the parameters mentioned in Sec. 
6.2 are shown in Fig. 6.1(a-f) in the form of probability density contoxirs superimposed on 
the potential energy contours in (i?, r) plane for a fixed value of 7 = 0.145 rad. Fig. 6.1(a) 
corresponds to the initial location of the WP. In about 8.08 fe the WP spreads along the 
asy mm p.trir stretch direction and in another 8.08 fe it develops structures. An analogous 
feature was found with the WP for the collinear geometry (c.f. Fig. (3.1)). The WP shows 
‘ a heavy build up of probability density along the (HeH"^ + H) channel (Fig.6.1(b)) and 
the subsequent time evolution (Fig. 6.1(c)) indicates high vibrational excitation of HeH 
and also of Hj. It can be seen from the nodal pattern of this time evolved WP that HeH 
is excited to the t/ = 5 and H^ to the t; = 8 quantum state. In 24.24 fe the WP gets into 



Figure 6.1: Probability density contours of the 
each box) showing the evolution in spac( 
{^*\^) of the WP retained on the grid at differei 


0.20 in (a), (b), (c), (d), (e) and (f) respectively. 
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Figure 6.3: Power spectrum in the energy range O.M.O eV computed by Fourier 
forming C(t) of the GWP. 


trans- 
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the well of the PES (Fig. 6.1(d)) and it reveals a nodal progression along the asymmetric 
stretch mode of the HeHa"- complex at later times as can be seen in Fig. 6.1(e) and 6.1(f). 
Since the masking fnnction (Eq. (3.13) )is activated at the end of both the channels we do 
not see any reflection of the WP from grid edges. 

The decay of the autocorrelation function (up to 2 ps) is shown in Fig. 6.2. The 
recurrence pattern is much simpler when compared to the colhnear results. The number 
of recurrences is much less and the amplitude dies off to zero within 1.35 ps. Therefore, 
it can be concluded that the dynamics of the system is relatively simpler in 3D as was 
anticipated by Sathyamurthy et al [255] from a study of the nature of the vibrationally 
adiabatic potential curves for the non-collinear geometries. 

The transition state spectrum obtained by Fourier transforming the autocorrelation 
function is shown in Fig. 6.3. Since the initial WP was centered on a high energy contour 
of the PES the peaks in the spectrum correspond to the high energy resonances. Fig. 
6.3(a) corresponds to the spectrmn obtained without using the Hanning window ftmction 
and 6.4(b) corresponds to the same obtained by using the latter. As can be seen from 
Fig. 6.3(b) the broad features of Fig. 6.3(a) are resolved into sharp features by the use 
of the window function. The spacing between the peaks in spectrum 6.3(b) are almost 
uniform indicating that the corresponding energy levels follow a regular progression. The 
computation of the resonances at other energies by varying the location of the WP and on 
a finer grid and a systematic analysis are being pursued. 

6.4 Summary and Conclusion 

The SQM has been used for investigating resonances in 3D (He,Hj) collisions. We have 
used the FFT method to evaluate the radial derivatives in the Hamiltonian and the Gauss- 
Legendre DVR for evaluating the angular derivatives. We have obtained preliminary re- 
sults on a coarse grid for the 3D He+Ht( J = 0) collisions which reveal the existence of 
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resonances. The dynamics of the system appears simpler in three dimensions. A more 
detailed analysis on a finer grid is presently in progress. 



Chapter 7 

Summary and Conclusion 


Time-dependent quantum mechanical wave packet methods have been used extensively 
in recent years in investigating problems relating to the area of TS spectroscopy and reac- 
tive scattering processes. The TDQM approach provides classical insight into the molecular 
encounter despite its quantum mechanical origin. Since a real chemical reaction occurs in 
the time domain this approach offers a direct link with experiment. The development of 
newer algorithms and faster computing machines has made it possible to study a wide 
variety of problems in the TDQM frame work. One major advantage of this theory is that 
a single WP calculation enables the prediction of experimental observables over a range of 
energy values. 

In chapter 1 of the thesis, we have presented an overview of the chronological develop- 
ment of the TDQM methodology. Its application to a variety of problems with a special 
emphasis on TS spectroscopy and reactive scattering is reviewed. The advantage of the 
TDQM method over the classical, semi-classical and TIQM methods is discussed. 

In chapter 2, the general grid method of solving the time-dependent Schrodinger equa- 
tion is presented. The representation of the initial WP on a grid and its evolution in space 
and time are described at length. We present the general collocation or pseudospectral 
method and its special case when orthogonal basis functions are used, in the context of 
evaluating This orthogonal collocation scheme leads to the Fourier method when 
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plane wave basis functions are used. The use of the fast Fourier transform (FFT) method 
in order to carry out the operation has been illustrated. We have also discussed the 
finite difference (FD) method in this regard and illustrated the superiority of the FFT 
method over the FD method. The use of the fast Hankel transform (FHT) for evaluation 
of the action of the radial part of the Laplacian operator in polar coordinates on is 
also discussed. Utilising the idea of the orthogonal collocation method, can also be 
evaluated locally when a basis set of some classical orthogonal polynomials is used. This 
results in a set of Gaussian quadrature sampling points in the discrete Hilbert space. A 
unitary transformation matrix is obtained in this situation which allows one to switch back 
and forth between the discrete local and non-local spaces. This method is known as the 
discrete variable representation (DVR) in quantum mechanics. This method is particularly 
useful in evaluating the action of the angular kinetic energy operator on ^ where the usual 
FFT technique ends up with singularity. This DVR method is discussed at length in this 
chapter. 

We have also discussed the various numerical time evolution schemes, viz., the SOD 
scheme, the SO method, the CP expansion scheme and the SIL method. We have paid spe- 
cial attention to the applicability of these schemes when the negative imaginary potential 
(NIP) is used to damp out the WP components near the grid edges. Our analysis shows 
that in presence of a NIP, the SOD scheme, is unstable,-the SO scheme results in an expo- 
nential damping of the WP components near the grid edges, the CP scheme requires the 
Hamiltonian to be renormalised in order to account for the shift of its eigenvalue spectrum 
in the complex energy plane and the SIL scheme also results in an exponential damping of 
the WP components and that the additional complications of the basis vector becoming 
non-orthogonal can be handled by the Gram-Schmidt orthogonalization scheme. 

We have presented the Fourier grid Hamiltonian (FGH) method for computing the 
bound state eigenvalues and eigenfunctions by utilising the variational method in the 



145 


TIQM frame work. The spectral quantization method for computing the bound states as 
well as the resonance eigenvalues by computing the power spectrum through the Fourier 
transform of the temporal autocorrelation function of the WP and the eigenfunctions is 
also described m this chapter. The numerical methods of computation of the vibrational 
state selected and energy resolved reaction probabilities through a time-energy mapping 
of the reactive flux of the WP and also the state-to-state reaction probabilities through a 
time-energy projection are presented. The strategy of computing state-to-state reaction 
probabilities by computing the S-matrix elements in the Mpller operator formalism is also 
presented. 

In chapter 3, we have presented the TS resonances and the dynamics of the collinear 
He + ^ HeH"^ ■(" H (HI) 

reaction and its isotopic variants when one of the two H atoms is replaced by a D on 
the MTJS PES. We have computed the TS spectrum of all three transients, HeHj , 
HeHD'*" and HeDH'*' by the spectral quantization method and the reaction probabil- 
ities {Ff{E)) for (HejHD"^) and (HejDH"^) collisions through a time-energy mapping of the 
reactive flux of the WP. Interestingly, we have found that the spectra for both HeHj and 
HeDH'^are dominated by a large number of TS resonances compared to that for HeHD'*', 
for which the threshold resonances seem to dominate. We have computed the eigenfunc- 
tions of the resonances and chairacterized them based on the local and hyperspherical mode 
description and assigned them appropriate quantum numbers. The analogy of the quantal 
resonances with the classical resonance periodic orbits is shown for the HeHf resonances. 
The life-times of some of the resonances are also reported. The TS spectrum of HeDH'*' is 
found to be densely packed when compared to that of HeHj and HeHD"^ and the corre- 
sponding dynamics is the slowest. The P^{E) values for (He,HD-") show a characteristic 
staircase-like structure which is related to the threshold resonances and those for (He,DH ) 
are highly oscillatory, in keeping with a densely packed TS spectrum. The markedly dif- 
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ferent dynamical behavior of the three systems emphasizes the role of kinematic effects 
in the system and particularly the complicated dynamics of coUinear (He,DH+) colUsions 
suggests the possible existence of "chaos" for which the classical trajectory dynamics would 
be worth investigating. The results obtained for the reaction (Rl) eire compared with the 
available TIQM results. 

In chapter 4, the TS resonances and the dynamics of the coUinear 
H" + H 2 -+ H 2 + H~ 

reaction on the SM PES have been reported. This system has shown the existence of a large 
number of TS resonances when compared to its neutral analog for which the resonances 
have been shown to be of threshold origin. Although the SM PES is very similar to 
the LSTH or the DMBE PES for the neutral analog the two systems differ significantly in 
terms of their dynamical behavior. The lf{E) values computed for reaction (R2) through 
the time-energy mapping of the reactive flux of the WP are highly oscillatory. We have 
explained the origin of some of the resonances therein on the basis of the quasibound states 
supported by the vibrationally adiabatic potentials. The eigenfunctions of the resonances 
reveal a matted structure at low energies and we have assigned a set of quantum numbers 
representative of the nodes along the hyperspherical radius (p) and angle (^). At higher 
energies the nodal progression along p disappears and results in the hyperspherical modes 
(progression along ) only. We are presently analysing the classical trajectoiy dynamics 
in order to explain the origin of the resonances in the system. 

The existence of a large number of narrow TS resonances in colUnear (He,Ht) coUisions 
and its known chaotic behavior in the classical domain inspired us to extract the signature 
of classical chaos in quantum domain. Although the term "chaos" seems to be well imder- 
stood in classical mechanics, its meaning and implication in quantum mechanics continues 
to be the subject of considerable debate. We have examined the short- range and long- 
range correlations in the eigenvalue spectrum of the TS of the system through a study of 
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the P( 5 ), and the AsCX), respectively. The spacing distribution statistics show an interme- 
diate behavior between the Poisson (regular) and the Wigner (irregular) limits. A similar 
conclusion is drawn from the A3(Ir)statistics, which show an intermediate distribution 
between the Poisson and the Gaussian orthogonal ensemble (GOE) limit. The behavior in 
the time domain has been examined by computing the ensemble averaged survival proba- 
bility, << by Fourier transforming the spectral autocorrelation function followed 

by averaging over the initial states and the Hamiltonian. << P(t) values when plotted 
against time t show the existence of the "correlation hole"— a typical quantum chaotic 
behavior. A quantitative comparison of « P{t) » with that obtained from the random 
matrix theory enables us to conclude that 70 % of the phase space exhibits chaotic behavior 
in collinear (He, Hj ) collision dynamics. We have also analysed the dynamical evolution 
of a coherent GWP located initially in different regions of phase space and computed the 
survival probability, power spectrum and the volume of phase space over which the WP 
spreads to illustrate the classical-quantal correspondence. The GWP initially centered at 
a chaotic point of classical phase space resulted in a spreading of the WP over the entire 
phase space. This leads to a large increase in the statistical entropy during the evolution 
of the WP. The survival probability values become totally uncorrelated in this chaotic 
situation. This has been docmnented in chapter 5 of the thesis. 

In chapter 6, we report on a preliminary investigation of the TS resonances in 3 D 
He+H^( J = 0 ) collisions on the MTJS PES. We have extended the SQM in three dimen- 
sions to compute the TS spectrum by Fourier transforming the temporal autocorrelation 
function, C'(f), of the initial WP written in terms of bond distances and bond angle. The 
action of the radial kinetic energy operators are evaluated by the FFT method and the ac- 
tion of the angular kinetic energy operator is evaluated by a 16 -point Gauss-Legendre DVR. 
We have obtained a few high-energy TS resonances of the systems on a (128 x 128 x 16 ) 
grid in (i?, r, 7) coordinates. The dynamics of the system in three dimensions looks simpler 
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tho-n that in the collinear geometry. The computation of other resonances by varying the 
location of the initial AVP in the interaction region of the PES and their systematic analysis 
on a finer grid are presently in progress. 
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Appendix 1 

I. Function of an Operator 

The eigenvalue equation of is given by; 

= Ek(l>k ( 1 ) 

If the basis functions [<j>k) form a complete orthonormal set then one can write resolution of 
the identity as 

1 Ipit >< hi (2) 

k 

and then the spectral resolution will be given by 

H = Y,Ek\h><<l>k\ ( 3 ) 

k 

Therefore any function /(^ can be expressed as 

'w 

k 

Hence the evolution operator will be given by 

(3) 

k 

II. The Conservation Laws 
II. 1. Norm Conservation 

The conservation of the norm of the wave packet during time evolution is a direct consequence 
of the uDitarity of the time eTOlutioo operator, OTT t = t/ 1 fr = 1. which has the additional 
requirements that E to be Hermitian. This can be illustrated as Mows: 

(®(()|«(f)) = (&«(0)|&®(0)) 

= ($(0)|fi'tp|®(0)> 

= («{0)|«(0)) 


( 6 ) 



In presence of NIPs H becomes complex and loses its Hermitian property. The unitarity of 
the evolution operator no longer remains valid and hence there will be no norm conservation. 
What is more important in this situation is, the efficient damping of the WP which reaches 
the grid boundary: 

The TDSE and its complex conjugate in presence of a NIP [— iVo] can be written as 


.^d<S 


2m 


+ y- iVo 




and 


ih 


dt 


*2 

_ v^+y+iVo 

2m 




Left multiplying Eq. (7) by and Eq. (8) by $ and then subtracting one arrives at 


a < > 2Vo 


dt 


h 


< > 


(7) 


( 8 ) 


(9) 


Solution of which results into 


< >t= >0 


( 10 ) 


an exponential damping of the norm of the WP. 

II.2. Energy Conservation 

Por a normalised WP the average energy at time t is given by 

{E)t = 

= {^{o)\u^mj\^o)) ( 11 ) 

H mw Hamiltoman commutes with the evolution operator then EU =VE, introducing this 
and with the help of the unitary property of t one can finaUy arrive at 

(E), = («(0)1^«(0)) 
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71115 rsm&iiis v&lid for sny form of fclio H flini ltoni&ri as long as it is time- indspondsnt and 
commutes with the time evolution operator. For time dependent Hamiltonians the evolution 
operator assumes a complex form and the commutation criterion no longer holds. 

II. 3. Conservation of the Time Reversal Symmetry 

Classical equation of motions are known to be invariant with respect to time reversal, which 
means all physical processes goes backward on interchanging the initial and final states 
of a system and its energy remains conserved. TDSE is a first order differential equation 
in time and is a primary physical law. As can be seen below TDSE is also symmetric with 
respect to time reversal. 

The TDSE and its complex conjugate in its differential form can be written as 


ih— = 

(13) 

— ih—;:— = 
dt 

(14) 


Now, let ^ = $(i2, t) is a solution of the first corresponding to the initial condition ^ 
for t = 0, then by substituting t = -f in the second one can also have = '^*{R,-t) 

corresponding to initial condition This symmetry is often called as the principle of 

microscopic reversibility or detailed balance. The time reversal operator i acting on the state 
function of a system results into its complex conjugate 

i(j>k = <f>k 


i is antilinear and antiunitary satisfies the equations 


ic(j)k = c*i<j)k 

(16) 


(IT) 


This symmetry operation leaves Hamiltonian unchanged and this symmetry exists as long as 

the Hamiltonian is real and Hermitian. 
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Appendix 2 

Rotational and Vibrational Excitation of HCl in Small 
Angle Scattering With He+ at 
Slab = 20-60 eV 

I. Introduction 

It has been demonstrated in the last several years that proton energy loss spectroscopy 
could be used to probe the vib-rotational inelasticity of target molecules (see for example, 
ref. 1,2). Recently it has been shown that highly energy resolved He^ could be used as a 
projectile in place of in such studies^. In particular, our recent paper^ illustrated the nicely 
resolved structures in the energy loss spectra of He"*" scattered in the forward direction due 
to collisions with HF and also how those structures could be explained in terms of rotational 
excitation in HF involving a change in rotational quantum number up to AJ = 35 with 
no indication of vibrational excitation. The observed propensity for large AJ transitions 
decreased with increasing collision energy, but increased with incretise in scattering angle. 
It was also shown that three dimensional(3D) quasiclassical trajectory (QCT) calculations^, 
treating HF as a rigid rotor and using a potential-energy surface (PES) based on the isotropic 
He‘''-Ne interaction at short range and ion-dipole (quadrupole) interaction in the long range 
reproduced all the essential features of the experimental results. In contrast the energy loss 
spectra of forward scattered He+ in collision with HCl reported in the present paper for Eiab in 
the range 20-60 eV reveal less structures, presumably due to the closer spaced rotational and 
vibrational levels for HCl (than for HF). While there are considerable large AJ transitions, 
there seem to be accompanying vibrational excitations as well, as was observed by Udseth et 
al’^ for H'^-HCl. In the spirit of our earlier paper^, we have modelled the He'^-HCl interaction 
at short range based on the isoelectronic He'-Ar® and in the long range using ion-dipole 
(quadrupole) interaction, with the damping functions introduced by Tang and Toennies’’ to 
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connect the two. We find that treating HCl as a rigid rotor is clearly inadequate to explain the 

observed expemnental results and it becomes necessary to include the vibrational coordinate 
of HCl in the model. 

The experimental details and results are presented in section n and the 3D QCT method- 
s^id results in section III. Su mm ary and conclusion are given in section IV. 

II. Experiment 

The experimental results were obtained using a crossed beam ion— molecule scattering ap- 
paratus described in detail elsewhere^’*. Usage of both electrostatic energy selection and 
analysis of the ions made it possible to achieve an energy resolution of AE (FWHM) = 30 
meV for a He'*' beam of = 20 eV, for example. To reduce kinematic broadening of 
the observed spectra the molecular beam was produced by expanding pure HCl gas at room 
temperature through a nozzle of diameter d = 70 fim with a stagnation pressure of po = 
190 mbar. The latter was adjusted to an attenuation of the ion beam of about 30% to avoid 
multiple collisions. Passing the beam through a skimmer located at 10 mm from the nozzle 
reduced the angular divergence to less than ±2°. Thus the overall spread in energy was kept 
at AE = 35 meV. This made it possible to resolve the final rotational states (J/ > 13) of 
HCl in its ground vibrational state (v = 0). The analyser-detector arrangement could be 
rotated around the axis of the secondary beam for measuring the differential cross section. 
The angular resolution in the laboratory firame was estimated at ±1°. 

Assuming thermal equilibrium, rotational relaxation during the expansion would result 
in a rotational temperature Trot « 90 K for the HCl molecules in the beam. The Boltzmann 
distribution at that temperature suggests that about 80% of the molecules would be in rota- 
tional states Jj < 2, corresponding to a mean initial rotational energy of about E(Ji) « 7.3 
meV. Due to the substantially large rotational constant of HCl (Be = 1-31 meV), however, 
the higher rotational levels may not relax fast enough to reach thermal equiUbrium. This 
would result in an overpopulation of higher states and thus a broadening of the structures 
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observed in the experimental spectra. 

Energy loss measurements were carried out at ©^,5 = 0" for five different He+ beam 
energies. Ejnj = 20, 25, 30, 40 and 60 eV. These energies correspond to Ecm = 18, -, 
27, -, and 54 eV respectively. Additional measurements were made at ©i^,^ = 2° and 3° 
at Eiai, =20 and 30 eV respectively. 

It is clear from the energy loss spectrum obtained at different Eiah for Oiab = 0°, shown 
in Fig.l, that there is substantial vib-rotational energy transfer resulting in AE of as ninrh 
as 2.4 eV. For v=0, rotational excitation up to J=27 could be identified. It can be seen from 
the figure that there is a dramatic decrease (by six orders of magnitude) in intensity of the 
energy loss spectrum, with increase in AE. In addition, there seem to be distinct features 
for large AE? for the different Eiab- 

Results of the energy loss spectra obtained for Ejob = 20 and 30 eV for ©joi = 2° and 3" 
respectively are included in Fig.2a and 2b, in that order. Perhaps the most important result 
obtained therein is that the intensity of the energy loss spectrum decreases with increase 
in Oiab , over the entire AB range investigated. 

III. Theory 

III.l. Potential-Energy Surface 

The long-range interaction of He'*' with HCl was obtained by the well known multipole mo- 
ment expansion upto fourth order^ : 

. mWPi(cos7) , Q(r)P2[cos'Y) , ^ 

V(^,r,7) = ^5 + ^3 +l 2R^ 2R^ ^ 

where R is the center-of-mass separation between He+ and HCl, r the intemuclear distance 

of HCl and 7 the angle between the HCl bond axis and R with 7 = 0 representing the 

approach of He+ at the Cl-end of HCl. Po and P2 are the Legendre polynomials, /x and 

Q are the dipole and quadrupole moments of HCl while Oo and are its isotropic and 
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Figure 2: Energy loss spectra for He‘^-HCl(v=0) at (a) Eiaj — 20 eV, ©jot 0 and 2 (b) 
Ejai =30eV, ©j^ = 0° and 3°. 




anisotropic polarizability values. Contributions due to the polarizability of He+ are negligible 
in comparison and hence have not been included in the present study. 

The ab initio calculated data for the variation of fx with r are available in the literature^® 
and they can be approximated by a linear relation 

li{r) =ar-b ( 2 ) 

for 1.8 < r < 3.4 a.u. as illustrated in Fig.3a. Although the quadrupole moment of HCl is 
known^^ at re, its variation with r is not known. However, dependence of Q on r is known for 
(see Fig.Sb). Assuming that the variation of Q with r in the neighbourhood of re would 
be the same for HCl and HF, we write 

Q{r) = Q{re)+c{r-r^) + d{r-r^f (3) 

The ab initio values of both parallel (a||) and perpendicular (a_) components of polariz- 
ability were reported by Gianturco and Guidotti^® for HCl at different intemuclear distances. 

The variation of a|| and (Xj_ with r is quadratic in nature. Since ao and 03 are related to 

a i and ctx through relations 


“0 = ^[2ai|H-ai] 

0:2 = 

we have assumed a similar type of variation with r for ao and ai (see Fig.3c).Thus 


(4) 


ao = Oof' ~ boT -h co 

ao = Oof' — h 2 T -f C2 (^) 


All the parameters a, b, Q(r,), c, d, o„. to, C, h, ej used in the study ate listed in Table 1. 
The short-range Het- - HCl interaction -was taken to be the same as that of the isoelectromc 
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Figure 3: (a) Variation of the dipole moment, (in D) of HCl with r-r^ (in a.u.). The points (*) 
indicate the ah initio results of ref.lO and the solid line is the straight line (eq.2) fit to those points for 
1-8 < r < 3.4 a.u.. (b) Variation of quadrupole moment, Q (in a.u.) of HF with r-re(in a.u.) indicated 
by the points(*). A similar variation is assumed for HCl in the neighbourhood, of its equilibrium 
bond distance and shown by the solid line (eq.3). (c) Variation of isotropic (oo) and anisotropic 
( 02 ) polarizability (in a.u.) of HCl with r-re (in a.u.). The points (*) are the values of the same 
calculated from the ah initio data (ref.l3) for the parallel (ay) and perpendicular (aj.) components 
-r — iioino- ATI 4 The solid lines are the quadratic fit to those points (eq.5). 




He'*' - Ar potential® because of the known nearly spherical nature of HCl. The short-ran<^e 
and the long-range potentials were connected smoothly using the damping functions of Tan^ 
and Toennies^ Representing the intramolecular potential of HCl by a Morse function the 
total - HCl interaction is given as 


V(R.r,7) = D, [l - 4- + Qfr)P2(cos7) ^ 

ao{r) a2(r)P2(cos7)' 


+ 


2R^ 


2W‘ 


( 6 ) 


with 






k=0 


kl 


(7) 


The parameters A and {3 and the Morse parameters De, and rg for HCl are included in 
Table 1. 

Potential-energy contours for the He'*' - HCl interaction for r=re of HCl are illustrated in 
Fig.4. There is a substantial well of depth 0.83 eV at R = 3.86 a.u. and 7 = 75°, Although 
there are no ab initio data available for He'*' - HCl interaction we can compare the above with 
a value of 0.863 eV obtained for He"*" - HF at R = 3.72 a.u. and 7 = 54°. These results are 
consistent with the known ab initio results^'* available for H"** - HF which reveal the presence 
of a deeper well at a shorter R revealing the greater penetrating power of H"^ than that of 
He+. 


III.2. The Rigid Rotor Model 

Since we could account for much of the energy transfer in He+-HF collisions treating HF as 
a rigid rotor (RR), we examined first the adequacy of the RR model for HCl-He-^- collision 

also. 

The trajectory code for the rigid rotor HCl - He+ collision was based on the algorithm 
outlined by Pattengill®. Fourth order Runge-Kutta method with a stepsize of 0.5387 fe was 



a 

-4.9504 eV a.u. 


6 

-0.1279 eV a.u.2 


Q(^e) 

75.8 eV a.u.3 

[11] 

c 

45.5392 eV a.u.^ 


d 

10.6127 eV a.u. 

- 

Oq 

54.424 eV a.u.^ 

- 

^0 

145.5842 eV a.u.^ 

- 

Co 

505.32684 eV a.u.^ 

- 

CL2 

44.8998 eV a.u.^ 

- 

^2 

93.0650 eV a.u.^ 

■ - 

C2 

-2.6123 eV a.u.'^ 

- 

A 

462.121 eV 

[6] 

/5 

1.7025 a.u.-^ 1 

[6] 

De 

4.6188 eV 

[17] 

Pm 

0.9886 a.u.“^ 

[17] 

Te 

2.4094 a.u. 

[17] 


Table 1. Potential parameters for He"^ - HCl. Units were 
adjusted to yield energy in eV. 
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Figiire 4: Potential - energy contours (in eV) for He'*' interaction with rigid HCl in its equilib- 
rium bond distance. For reference the size of H and Cl are indicated, based on their covalent 
radii. Distances are in a.u.. 


used to solve the Hamilton’s equations. The HCl molecule was taken to be initially in its 
ground rotational state (Ji = 0). Trajectories were started at R = 35 a.u. where the He+ - 
HCl interaction was negUgible. The initial value of the center-of-mass relative translational 
energy (Etrans) was fixed at 18, 27, or 54 eV corresponding to Biab— 20, 30 and 60 eV 
respectively. The initial impact parameter (b) was sampled uniformly as follows; 

3.5 (0.1) 9.0 (0.5) 29.0 a.u. at 18 eV; 

3.5 (0.1) 10.0 (0.5) 26.5 a.u. at 27 eV; 

3.4 (0.2) 10.0 (1.0) 23.0 a.u. at Etrans= 54 eV. 

All other variables were sampled randomly. Typically, 2000 trajectories were rim at each 
value of b. The minimum and the maximum values of b were chosen by examining the final 
rotational state (J/) of HCl and the scattering angle 

In view of the experiments having been carried out for ©ia 5 = 0° with a resolution of 
~ 1°, which corresponds approximately to ©c.m.= 0°-l°, we were interested mainly in the 
trajectories leading to 0c.m.= 0° - 1°- The final rotational state distributions were arrived at 
by using the standard histogram procedure®. However, it must be added that for b > 26.5 
a.u. at Efrans= 27 eV, for example, some of the trajectories led to 0.5 <. Jy < 1.0 indicating 
that the 0 — > 1 rotational transition was classically forbidden^®. Therefore the collisions were 
assumed to be elastic at b > 26.5 a.u. Similarly the Ji = 0 J/ = 2 transition was found to 

be classically forbidden at b > 20.5 a.u. and trajectories with 1.5 < J/ < 2.0 were included 

under J/ = 1-0. A similar procedure was adopted for other Jy values also. 

The state-to-state inelastic transition probabilities P were computed as follows: 


where is the number of trajectories resulting in Jy ± 0.5 (vide supra) and Ntot the 

total number of trajectories computed at a particular b. The dependence of on b for 

diflFerent Jy values is shown in Fig.5, for Etrans= 27 eV. As has been observed for other systems 
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Figure 5: Opacity fiinctioiis for different J/, for rigid rotor HCl(Ji=0)-He'*' collisions at Etrans 
— 27 eV for scattering angle (©c.m.) ia the range 0° - 1°. The orbital angular momentum (1) 
in units of ft are shown in the upper abscissa of Fig. (a) corresponding to the b values in the 
lower abscissa. 





( for example, see ref.4 ), small AJ transitions occur at large b while large AJ transitions 
occur at small b. It is worth emphasizing that the amount of orbital angular momentum ( 
Icri = /X Vrei b, where /i is the reduced mass of He-*- - HCl and Vrd the relative velocity of 
the collision partners ) is two orders of magnitude leirger ( see Fig.5 ) than the amount of 
rotational angular momentum transferred. 

The state-to-state differential cross sections for 0c.m.== - 1° were computed using the 

relation: 


^ (©c.m.) 


SZT10C.771. 


{b, 0c.m.) 


b Ah 

A0c.Tn. 


and the results are plotted as a function of the amount of energy transferred ( AE) for J, 
= 0 at different Etrans inFig.6. The largest energy transferred is about 1500 meV, for Etroni= 
IS eV, for which the experimental results indicate an energy loss upto 2500 meV. In principle, 
there could be trajectories resulting in energy loss greater than 1500 meV but they would 
originate at small b and sample regions in (R,7) space for which the interaction potential used 
is of questionable validity. Hence we have not computed such trajectories. Also, our 3D QCT 
results have already given an intensity variation of five orders of magnitude. Reproduction 
of larger AE results would require an additional order of magnitude increase in the number 
of trajectories computed. 

A comparison of the simulated energy loss spectra for Etnins= 18, 27 and 54 eV with 
the corresponding experimental results shown in Fig.7 reveals that the agreement between 


theory and experiment is excellent for AE up to 1500 meV at the lowest Etrons used. The 
range of energy loss for which the agreement between the predicted and the observed spectra 
is satisfactory, decreases dramatically with increase in Etrans- 

The maximum value of observed in our trajectory studies also decreases with increas- 
ing Etrans as can be seen from Fig.6 : 34, 27 and 18 at 18, 27, and 54 eV respectively. This 
is definitely due to a decrease in the interaction time with increase in Et^on* and not to the 
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Figure 6: Differential cross sections for scattering in the range Scm. - 0“ - 1°, for rigid rotor 
HCl(Ji=0)-He''' collisions at Etrans = 18, 27 and 54 eV. ■ 
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Figvire 7: Simulated energy loss spectra compared with the experimental results for rigid 
rotor HCl(Ji=0)-He+ collision for 0cm. = 0° - 1° at Etron. = 18, 27 and 54 eV. 
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lack of availability of angular momentum in the collisions. Similar results were obtained for 
He - HF collisions and they were in agreement with experiment. 

The dramatically large discrepancy between theory and experiment for HCl-He+ at large 
Etrans could be due to several factors: 

(i) inadequate representation of He+ - HCl interaction 

(ii) use of classical trajectory approach for investigating the dynamics 

(iii) neglection of contributions from nonzero Jj states of HCl 

(iv) poor energy resolution in He'*' beam and 

(v) inadequacy of the rigid rotor model for HCl. 

The fact that a classical mechanical treatment using a similar potential ( short range 
mimicking He”^ - Ne and asymptotically accurate long range ) led to reasonable agreement 
between theory and experiment for He+ - HF suggests that the first two factors are unlikely 
reasons for the discrepancy. Classical trajectory calculations for Jj = 1 and 2 of HCl led to 
very similsir inelastic cross sections as for Ji = 0 at Etrans= 18 eV as illustrated in Fig.8 and 
this was the case for He''' - HF as well. The known energy spread (35 - 65 meV ) in the He'*' 
beam is too httle to account for a discrepancy in AE of the order of 1000 meV. Therefore, 
it is reasonable to conclude that the rigid rotor model for HCl is not adequate for explaining 
the large energy transfers at high Etrans- 

Also it may be pointed out that both the rotational and vibrational constants axe smaller 
for HCl than for HF. This would suggest a larger vibrational excitation cross section in the 
case of HCl as a target as was found by Udseth et al^ at lower resolution. This also means 
that with increasing Etrans the vibrational excitation cross section would become larger. 
This is indeed what is seen in our experiments, cf. the peak at v = 1. The observed peaks 
in the energy loss spectra coincide with the calculated positions for vibrational transitions of 

HCl accompanied by rotational energy transfer. 

Still there remain some structmes in the measured spectra, which cannot be expired 




Figure 8: Differeutial cross sections for scattering in the range ©am. — 0° - 1 , for rigid rotor 
HCl(Ji=0,l,2)-He+ collisions at = 18 eV. 
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fully at the present level of investigation. These may be due to the excitation of rather high 
(J/ 15) rotational states in a given vibrational state v which involve the sa_me amount of 

energy transfer as the next higher vibrational transition and cannot bcrcscolved. There may 
also be some combination of v and J f which are favourably excited in foiwsrJl scattering. 
These structures do not appear in spectra recorded under larger (> I'*) sca.ttesiing angles. 
Electronic excitation and charge transfer cannot account for these features au ithey involve 
energy transfers, greater than 5 eV, not covered by our experiments. Also it may be added 
that the differences arising from the different rotational and vibratonaltmraosifcions for the 
different isotopomers of HCl occuring in nature (^®Cl : ^"^Cl = 3:1) would btarcdly’ account for 
differences of the order of a few meV, far below the resolution of the lecorcdedl sijectra. 

Therefore it becomes clear that the RR model is inadequate for HCl ia exq)laining the 
experimental results and it becomes necessary to include the vibrational cooidnnate of HCl 
in the simulation. 

III.3. T?he Vibrating Rotor Model 

In the vibrating rotor model the initial intemuclear distance of HCl was select ed randomly 
using the relation’^® 


r = Te 


_ ± In r {-2A) [B + VB'-4ACsiji (2r4)] 

PM ^ 



( 10 ) 


where are the Morse parameters defined above, ^ is a randotm muHuber varying 

between 0 and 1 and 


The other variables were selected randomly using the standard pxoce(Lure[5)]. The fourth 
order Runge-Kutta method was used to solve the Hamilton’s equaticsns off motion with a 
step-size of 0.10775 fs. Since the long range part of the potentiad plays sa *oMDinant role in 
forward scattering we had to start our trajectory at a large R value «f 50 a. u., where the 


interaction potential is practically negligible. 

The HCl molecule was initially taken in the Vi = 0, Ji 


state- Th^e Lnitial relative 
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translational energy was taken to be the same as in the rigid rotor case. The impact parameter 
was sampled uniformly as follows; 

3.4 (0.1) 4.0 (0.2) 29.0 a.u. at Etran. = 18 eV; 

3.4 (0.1) 4.0 (0.2) 29.0 a.u. at Etron, = 27 eV; 

3.4 (0.2) 4.0 (0.2) 23.0 a.u. at Etrans = 54 eV. 

and all the other parameters were sampled randomly as before. 

We have analysed the trajectories with the scattering angle &c.m. — 0-1° corresponding 
1-0 ^lab — 0° as well as ©c.m. == 2-4° corresponding to Qiaj, = 3°. The final rotational and 
vibrational quantum numbers were calculated following the standard histogram procedure^. 
The inelastic transition probabihties were calculated by: 


^Tat 


( 11 ) 


where is the number of trajectories originating from an initial (uj, Ji) state and 

leading to {vf ± 0.5, Jf ± 0.5). Ntu is the total number of trajectories computed at that 
particular impact parameter b. 

The state-to-state differential cross sections for 0c.m.= 0°-l“ computed using the 
relation: 


c.m.) 


Pvi,Ji-*Vf,Jf (5, 0C.771.) 


b Ab 

A0c.m. 


( 12 ) 


d uj- sinQc.m. 

The results plotted in Fig.9 show that with increase in Etrons the amount of energy transferred 
decreases. For Vf = 0, the extent of rotational energy transferred is approximately the same 
as obtained from the rigid rotor model. But there is significant energy transferred into Vf 


= 1 level. The overall energy transferred as obtained from QCT calculations is shown in 
the form of simulated spectra and compared with the experimental results in Fig.lO. The 
agreement is very good for Et„„5=18 eV, for 1600 meV. The agreement between theory 

and experiment worsens with increase in Etrona- While it is tempting to suggest that the 
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Figure 9: Differential cross sections for the vibrating rotoj. HCl(vi=0,Ji=0)-He+ collision hr 
three different Etransi (a) eV (b) 27 eV and (c) 54 eV, fQj. 0 ^^ = 0 °-l° (soHd lines) and 
2°-4° (dashed lines). The maximum rotational excitation diff^ent final vibrational 

levels of HCl(v/) are indicated by the arrows in the Fignjg 




Figtire 10: Simulated energy loss spectra compared with the experimental results for the 
vibrating rotor HCl(vi=0,Ji=0)-He+ coUision for 0cm. = 0°-!°. 



discrepancy is due to inadequate representation of the potential, we considered the possibiUty 
of the range of ©c.m.being larger than 1°. That is, we included trajectories with ©c.m.= 0-2’ 
in the simulation. The agreement improves marginally at all energies as shown in Fig.ll. We 
have also computed the energy loss spectra for collisions resulting in ©c.m.= 3 ± 1° and they 

are in very good agreement with the experimental results as illustrated in Fig. 12 for Ejran 5 = 
18 and 27 eV. 

The striking discrepancy that persists at higher Etrans in all these pictures is that the 
QCT calculations underestimate the extent of large AE transfer, arising, presumably due to 
large Au transitions. This could arise , from the fact that the PES used does not satisfactorily 
depict the r-dependence of the potential at short range. Because of large vibrational level 
spacings in HF, the rigid rotor model was adequate for HF and the He'''-Ne interaction served 
as a reasonable representation of the potential at short R, in the case of He'-BP. For He'^-HCl. 
on the other hand, the inclusion of the r-coordinate becomes essential, but, unfortunately, 
there is little information available, experimental or theoretical, for the vibrational couphng 
of the short-range interaction. 

rV. Summciry and Conclusion 

Energy loss spectra obtained for He^ collisions with HCl at ^lab— 20-60 eV for forward 
scattering of He'*' revealed a dramatic (six order of magmtude) decrease in intensity with 
increase in AE. The spectrum obtained at the lowest Etrans (18 eV) could be accounted 
for, almost quantitatively, for AE up to ~ 1500 meV using 3D QCT calculations treating 
HCl as a rigid rotor. By including the vibrational coordinate of HCl in the QCT calculation, 
the agreement improves, till ~ 1700 meV. By allowing a slightly larger (0-2 ) margin for the 
scattering angle, experimental results upto ~ 2000 meV could be reproduced by theory. 

With increase in Eiaj,,' there is very little variation in the extent of energy transferred, 
although there seems to be some propensity for large Av transfer. The agreement between 
theory and experiment is good for small to medium AE even at large Etr<m» • But 3D QCT 


Intensity (count 



-190 190 570 950 1330 1710 7090 2470 

c.m. Energy Loss (meV) 


Figure 11: Same as Fig.lO for 0c.m. 
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calculations fail to reproduce the large tlE at large Etrans , presumably due to inadequate 
representation of the r-dependence of the potential. Unfortunately there is no experimental 
or theoretical input available for this aspect of the He'''-HCl interaction. Hopefully, ab initio 
calculations might provide insight into the same, in the near future. 
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Addendum 


1. p. 90, second line from the bottom please include ref. 200. 

2. Ref. 137 should be read as; 

T. Seideman and W. H. Miller, J. Chem. Phys. 96, 4412 (1992). 

3. Ref. 280 should include 

Z. Top and M. Baer, Chem. Phys. 25, 1 (1977). 

4. p. 11, Below Eq. (2.14) G should be read as 

^in “ 9n{p^i) 


5. p. 20, Eq. (2.47) the summation should be read over n. 

6. p. 54, Eq. (3.13) is replaced by 


/(Xi) = sin 


^ i,^mask “I” ^Xjnask ^i) 


^^mask 


^ ^maak 


Xmaak is the point at which the masking function is initiated and AXmaaki= Xmax — 
Xmaak) is the width over which the function decays from 1 to 0. 

7. p. 81, line 5. More recent calculations revealed a slight difference in the bound 
state energies of HeHD'*’ and HeDH'''. The energy values are -0.1499, -0.0655 eV for 
the former and -0.1514, -0.0664 eV for the latter. 

8. p. 95, five line up from the bottom, "The upper potential" should be read as 
"The upper bound of the potential". 

9. p. 130, line 11, (b) and (c) should be read as (a) and (b). 



